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Abstract 



This study derives geometric, variational discretizations of continuum theories arising in fluid dy- 
namics, magnetohydrodynamics (MHD), and the dynamics of complex fluids. A central role in these 
discretizations is played by the geometric formulation of fluid dynamics, which views solutions to the 
governing equations for perfect fluid flow as geodesies on the group of volume-preserving diffeomorphisms 
i i of the fluid domain. Inspired by this framework, we construct a finite-dimensional approximation to the 

i-C diffeomorphism group and its Lie algebra, thereby permitting a variational temporal discretization of 

geodesies on the spatially discretized diffeomorphism group. The extension to MHD and complex fluid 
fH flow is then made through an appeal to the theory of Euler-Poincare systems with advection, which pro- 

vides a generalization of the variational formulation of ideal fluid flow to fluids with one or more advected 
parameters. Upon deriving a family of structured integrators for these systems, we test their performance 
via a numerical implementation of the update schemes on a cartesian grid. Among the hallmarks of these 
new numerical methods are exact preservation of momenta arising from symmetries, automatic satisfac- 
tion of solenoidal constraints on vector fields, good long-term energy behavior, robustness with respect 
to the spatial and temporal resolution of the discretization, and applicability to irregular meshes. 
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1 Introduction 

Fluids, magnetofluids, and continua with microstructure are among the most vivid examples of physical 
systems with elaborate dynamics and widespread appearance in nature. For systems of this type - whose 
motions are governed by collections of coupled partial differential equations that include the Navier-Stokes 
equations, Maxwell's equations, and continuum counterparts to Euler's rigid body equations - numerical 
algorithms for simulation play an indispensable role as predictors of natural phenomena. The aim of this 
paper is to design a family of structured integrators for magnetohdyrodynamics (MHD) and complex fluid 
flow simulations using the framework of variational mechanics and exterior calculus. 

Naive discretizations of physical theories can, in general, fail to respect the physical and geometric struc- 
ture of the system at hand [24] . A classic example of this phenomenon arises via the application of a standard 
forward-Euler integrator to a conservative mechanical system, where some of the most basic structures of 
the continuous flow - energy conservation, volume preservation, and the preservation of any quadratic or 
higher-order invariants of motion like angular momentum - are lost in the discretization. Even high-order in- 
tegration schemes, including popular time-adaptive Runge-Kutta schemes for ordinary differential equations, 
are prone to structure degradation, unless special care is taken to design the integrator with the appropriate 
goals in mind [36]. 

Conventional numerical schemes for MHD and complex fluid flow bear similar defects. It is well known to 
the numerical MHD community, for instance, that failure to preserve the divergence-freeness of the magnetic 
field during simulations can lead to unphysical fluid motions [13]. A vast array of MHD literature over the 
past few decades has been devoted to this issue. Proposed solutions include divergence-cleaning procedures, 
the use of staggered meshes, the use of the magnetic vector potential rather than the magnetic field, and 
even modification of the MHD equations of motion themselves [14, 32, 41, 45]. 

For a variety of physical systems, the problem of structure degradation in numerical integration may be 
remedied through the design of variational integrators [31]. Such integration algorithms made their debut 
in mechanics, where a discretization of the Lagrangian formulation of classical mechanics allows for the 
derivation of integrators that are symplectic, exhibit good energy behavior, and inherit a discrete version 
of Noether's theorem guarantees the exact preservation of momenta arising from symmetries [28, 36]. An 
extension of these ideas to the context of certain partial differential equations may be made through an 
appeal to their variational formulation [34] , and the associated integrators can often be elegantly formulated 
in the language of exterior calculus on discrete manifolds [2, 3, 9, 17, 25]. 

The fields of magnctohydrodynamics and complex fluid dynamics provide intriguing arenas for the design 
of structure-preserving integration algorithms. The former field lies at the confluence of two major domains of 
the physical sciences - fluid dynamics and electromagnetism - and its equations of motion comprise the union 
of two celebrated systems of equations in modern physics: the Navier-Stokes equations, describing the motion 
of fluids, and Maxwell's equations, describing the spatial and temporal dependence of the electromagnetic 
field [23]. The theory of complex fluid flow is equally interdisciplinary, entwining fluid dynamics with 
continuum versions of Euler's rigid body equations. 

Progress in the Structured Integrators Community. Conveniently, structured integrators for fluid, 
electromagnetism, and rigid body simulations have been subjects of active research as we briefly review. 

For fluids simulation, Perot et al. [40, 46] and Mullen et al. [38] have developed time-reversible integrators 
that preserve energy exactly for inviscid fluids. Elcott and co-authors [19] have proposed numerically stable 
integrators for fluids that respect Kelvin's circulation theorem. Cotter et al. [16] have provided a multisym- 
plectic formulation of fluid dynamics using a "back-to-labels" map, i.e., the inverse of the Lagrangian path 
map. More recently, Pavlov and co-authors [39] have derived variational Lie group integrators for fluid dy- 
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namics by constructing a finite-dimensional approximation to the volume-preserving diffeomorphism group 
- the intrinsic configuration space of the ideal, incompressible fluid. The latter advancement shall serve as 
the backbone of our approach to the design of integrators for MHD and complex fluid flow. 

In the electromagnetism community, Stern and co-authors [42] have derived multisymplectic variational 
integrators for solving Maxwell's equations, proving, in the process, that the popular Yee scheme [44] (and 
its extension to simplicial meshes [10]) for computational electrodynamics is symplectic. The integrators of 
Stern et al. [42] are formulated in the framework of discrete exterior calculus, giving the added benefit that 
the integrators may be easily made asynchronous. 

Structured rigid body integrators have a longer history, due in large part to the simpler nature of 
finite-dimensional systems. Among the many milestones in structured rigid body simulation are Moser 
& Veselov's [37] integrable discretization of the rigid body, Bobenko & Suris's [8] integrable discretization 
of the heavy top, and Bou-Rabee & Marsden's [12] development of Hamilton-Pontryagin-based Lie group 
integrators. For a more comprehensive survey, see reference [12]. 

In this report, we combine techniques from the structured fluid, structured electromagnetism, and struc- 
tured rigid body communities to design a family of variational integrators for ideal MHD, nematic liquid 
crystal flow, and microstretch continua. The integrators we derive exhibit all of the classic hallmarks of vari- 
ational integrators known to the discrete mechanics community: they are symplectic, exhibit good long-term 
energy behavior, and conserve momenta arising from symmetries exactly. Moreover, their formulation in the 
framework of discrete exterior calculus ensures that Stoke's theorem holds at the discrete level, leading to 
an automatic satisfaction of divergence-free constraints on the velocity and magnetic fields in the resulting 
numerical schemes. 

Layout. This paper consists of five main sections. In Section 2, we describe the geometric formulation 
of ideal fluid flow, and we explain the role played by the diffeomorphism group in this formulation; we 
summarize some key Lie group theoretic aspects of the diffeomorphism group, and proceed to construct a 
finite-dimensional approximation of the diffeomorphism group in the manner laid forth by Pavlov and co- 
authors [39]. In Section 3, we present the theory of Euler-Poincare systems with advected parameters, which 
provides the variational framework for all of the subsequent continuum theories presented in this paper. We 
then derive a variational temporal discretization of the Euler-Poincare equations with advected parameters. 
In Section 4, we state precisely the geometric formulation of ideal fluid flow and proceed to discretize it using 
the tools developed in Sections 2-3. We then discretize three continuum theories: magnetohydrodynamics in 
3D, as well as nematic liquid crystal flow and microstretch fluid flow in 2D. We present these discretizations 
as methodically and comprehensively as possible in order to highlight the systematic nature of our approach. 
In Section 5, we specialize to the case of a cartesian mesh and record the cartesian realizations of the 
numerical integrators derived in Section 4. Those readers most interested in computational implementation 
may wish to proceed directly to the end of this section for a concise catalogue of our novel numerical schemes. 
Finally, in Section 6, we implement our structured integrators on a variety of test cases adapted from the 
literature. We focus primarily on our MHD integrator since, relative to complex fluid dynamics, the field 
of computational MHD is replete with well-established numerical test cases and existing integrators for 
comparison. We show numerically that our integrators exhibit good long-term energy behavior, preserve 
certain conserved quantities exactly, respect topological properties of the magnetic field that are intrinsic to 
ideal magnetohydrodynamic flows, and are robust with respect to the spatial and temporal resolution of the 
grid. 

Our exposition is largely self-contained, but assumes a working knowledge of Lie groups and Lie algebras. 
For the reader's convenience, we give a brief summary of those aspects of Lie theory most relevant to our 
study in Appendix A.l. 
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2 The Diffeomorphism Group and its Discretization 

Pioneered by Arnold [4], the variational formulation of ideal fluid flow stems from the recognition that the 
governing equations 

dv 

— + (v-V)v = -Vp (2.1) 
V-v = (2.2) 

for ideal, incompressible fluid flow describe geodesic motion on the group Diff vo i(Af) of volume-preserving 
diffeomorphisms of the fluid domain M . Equivalently, in the language of mechanicians, (2.f-2.2) are Euler- 
Poincare equations on Diff vo i(M) with respect to a Lagrangian given by the fluid's kinetic energy 



|v|| 2 dx. (2.3) 

M 



This formulation has deep consequences in the analysis of fluid dynamics [33] and, as recently demonstrated 
by Pavlov and co-authors [39], can provide a powerful framework for numerical discretizations of fluid flows. 

For our purposes, the key pieces of insight deserving emphasis here are (1) that the configuration space 
of the ideal fluid is a Lie group, and (2) that the equations of motion (2.1-2.2) on this group are variational. 
The approach of this paper will be to construct a finite-dimensional approximation to the group Diff vo i(Af) 
in the manner of Pavlov and co-authors [39] and design variational Lie group integrators on the spatially 
discretized diffeomorphism group. An extension to MHD and complex fluid flow will later be made through 
an appeal to the theory of Euler-Poincare systems with advection, which provides a generalization of the 
variational formulation of ideal fluid flow to fluids with one or more advected parameters. 

We devote this section to a study of the geometry of the diffeomorphism group, followed by the construc- 
tion of a spatially discretized diffeomorphism group. 



2.1 The Continuous Diffeomorphism Group 

Let M be a smooth manifold, hereafter referred to as the fluid domain. The volume-preserving diffeomor- 
phism group Diff vo i(A/) consists of smooth, bijective maps if : M — > M with smooth inverses. The group 
multiplication in Diff vo i(M) is given by function composition. 

The Lie algebra of Diff vo i(M) is Xdiv(M), the space of divergence-free vector fields tangent to the bound- 
ary of M. Fixing a volume form dx on Af, the space dual to Xdw(M) may be identified with fi (M)/eJf2°(M), 
the space of one-forms on M modulo full differentials, under the pairing (•, •) : Xdiv(M)* x X<u v (M) — > K 
given by 




for any w b € J7 1 (M), v £ X<iw{M). Here, \vr] denotes the coset of one-forms in fi 1 (Af) / dfl° (Af) with repre- 
sentative w b . For easier reading, we will suppress the brackets when referring to cosets in fi 1 (Af) / dfl° (M) 
for the remainder of this paper. 

Adjoint and Coadjoint Actions. Let ip € Diff vo i(M), u,v e £ div (M), and w b e fl 1 (M) / dfl° (M) . The 
adjoint and coadjoint actions on Xdi V (Af) and its dual are, respectively, the pushforward of vector fields and 
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the pushforward of one-forms: 

Ad v v = (p*v (2.5) 
Ad;- lW b = <^V. (2.6) 

The infinitesimal adjoint and coadjoint actions on X^wiM) and its dual are given by Lie differentiation: 

ad u v = — £ u v (2.7) 
ad^w b = £ u w b . (2.8) 

Note the sign of (2.7); it says that the Lie algebra bracket on Xdi v (M) is minus the standard Jacobi-Lie 
bracket of vector fields. 



2.1.1 The Group Action on Scalar Fields 

The group Diff vo i(M) acts naturally from the right on J-{M) = f2°(M), the space of scalar fields (zero-forms) 
on M, via the pullback: 

f-<p = <p*f (2.9) 

for any / £ F(M), ip G Diff vol (Af). 

The induced infinitesimal action of an element u of Xdw{M) on J-(M) is given by Lie differentiation: 

/ • u = £ u f. (2.10) 

We close this subsection with two remarks regarding the nature of the action of Diff vo i (M) on scalar fields 
that shall motivate the definition of our discrete approximation to the volume-preserving diffeomorphism 
group. First, any cp e Diff vo i(M) clearly preserves constant functions / € J-(M): 

f><p = f if / = const. (2.11) 

In addition, a theorem attributed to Koopman [30] states that such diffeomorphisms preserve inner products 
of scalar functions /, h: 

{f-tp,h-<p)=(f,h) (2.12) 

Here, the inner product is taken to be the standard inner product of scalar fields. Note that the latter 
property relies on the fact that tp is volume-preserving, while the former property does not. 



2.2 The Discrete Diffeomorphism Group 

Given a mesh M on the fluid domain M with cells C,, i = 1, 2, . . . , N, define a diagonal N x N matrix f2 
consisting of cell volumes: Qu — vol(Cj). To discretize the group Diff vol (M), define 

V{M) ={q£ GL(N)+ \ = 1 Vi, q T VLq = ft}, (2.13) 

3 

the group of fi-orthogonal, signed stochastic matrices. Elements of X>(M) are referred to as discrete diffeo- 
morphisms. Later we shall see that matrices in X>(M) respect properties (2.11-2.12) in a discrete sense. 
The Lie algebra of 2?(M), denoted by Z)(M), is the space of fi-antisymmetric, row-null matrices: 



D(M) q{(N) | A a = ATn + ^A = 0}. 

3 



(2.14) 
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Elements of 0(M) are referred to as discrete vector fields. 

Remarkably, the space dual to 0(M) may be identified with a discrete analogue of Cl 1 (M) / dCl° (M) under 
a pairing given by a well-known matrix inner product. The following definitions and the subsequent theorem 
make this statement precise. 

Definition 2.1. (Discrete Zero-Forms.) A discrete zero-form (scalar field) is a column vector F € M. N . 
The components of such a vector F are regarded as the cell averages of a continuous scalar field f € F{M), 
i.e. Fi — J c .fdx J vol(Cj). The space of discrete zero-forms is denoted fl d (M). 

Definition 2.2. (Discrete One-Forms.) A discrete one-form is an antisymmetric matrix C b € so (AT). 
The space of discrete one-forms is denoted fi^(M). 

Definition 2.3. (Discrete Exterior Derivative of Zero-Forms.) The discrete exterior derivative is the 

map d : Q d (M) — > fl d (M) taking a discrete scalar field F to the discrete one-form dF whose entries are given 
by 

(dF) ij =F i -F j , i.j 1.2 V. (2.15) 

The image of d is denoted dO°(M) and its constituents are referred to as discrete gradients or full discrete 
differentials. 

Theorem 2.4. (The Space Dual to d(M).) The space dual to D(M) may be identified with Q, d (M)/dQ d {M), 
the space of discrete one-forms on M modulo full discrete differentials, under the pairing (•,•) : 5(M)* x 
D(M) — > K given by the ^.-weighted Frobenius inner product 

= Tr(C bT flB). (2.16) 

Here, B is a discrete vector field and [C b ] denotes the coset of discrete one-forms in Q d (M)/dQ d (M) with 
representative C G f2^(M). 

Proof, ft is well known that the standard unweighted Frobenius inner product of matrices permits the 
identification of so(N) with its dual. Since the set f20(M) := {VlB\B <G D(M)} constitutes the subspace of 
matrices in so (TV) that are row-null, it suffices to show that df^(M) is the orthogonal complement to ild (M) 
in so(iV) with respect to the standard unweighted Frobenius inner product of matrices. 

To show this, observe that f20(M) has codimension N— 1 inso(TV). Indeed, any D e f2t>(M) is constrained 
by a system of N — 1 independent equations: 

5^Dy=0, » = 1,2,...,JV-1. 

3 

(The nullity of the N th row of D follows from the N — 1 equalities above together with the antisymmetry of 
D.) 

On the other hand, the space dVl d {M) of discrete gradients has dimension N — 1, since any such quantity 
is defined by the N components F\, F 2 , . . . , Fn of a discrete zero- form F, modulo a constant additive factor. 
Moreover, the standard unweighted Frobenius inner product Tr((dF) T £>) of any D € f2t>(M) with a discrete 
gradient dF € dfl d (M.) is zero by the row-nullity and antisymmetry of D. Hence, the orthogonal complement 
to 05 (M) in so(.ZV) with respect to the standard unweighted Frobenius inner product coincides with the space 
of discrete gradients. ■ 

As before, we will suppress the brackets when referring to cosets in Q d (M)/dQ d (M) for the remainder of 
this paper. 
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Adjoint and Coadjoint Actions. Let q £ V(M), A,B £ 3(M), and C b £ 3(M)*. The adjoint and 
coadjoint actions on 0(M) and its dual are given by matrix conjugation: 

Ad q B = qBq- 1 (2.17) 
AdJ-iC* = q&Qq^n- 1 . (2.18) 

The infinitesimal adjoint and coadjoint actions on 0(M) and its dual are given by matrix commutation: 

ad A B=[A,B] (2.19) 
ad* A & = -[A, C^fi" 1 . (2.20) 

Formulae (2.17) and (2.19) are universal for matrix groups; formulae (2.18) and (2.20) are dictated by the 
pairing (2.16). 



2.2.1 The Group Action on Scalar Fields 

A right action of the group D(M) on f2°(M) can be constructed in the following manner: 

F-q = q- 1 F (2.21) 

for any F £ f2^(M), q £ T>(M). This action is merely the multiplication of a vector F by a matrix the 
inversion of q ensures that the action is a right action rather than a left action. 

The induced infinitesimal action of an element A of Z)(M) on fl3(M) is again given by multiplication: 

F ■ A — -AF. (2.22) 

With respect to the f2- weighted Euclidean inner product (F,H) :— F T VlH on Mr, the group action de- 
scribed above satisfies a pair of discrete analogues of properties (2.11-2.12) of the continuous diffeomorphism 
group action on scalar fields. Specifically, any q £ D(M) preserves constant functions: 

F-q = F if F = (c,c, ...,c) T G n° d {M),c £ R (2.23) 

In addition, such discrete diffeomorphisms preserve inner products of discrete scalar fields F, H: 

(F-q,H-q} = (F,H). (2.24) 

The latter property relies on the fact that q is f2-orthogonal, while the former property follows from the fact 
that q is (signed) stochastic. 

Properties (2.23-2.24) illuminate our choice (2.13) of a finite-dimensional approximation to the diffeo- 
morphism group. 



Correspondences. Table 1 summarizes the correspondences between the continuous and discrete diffeo- 
morphism groups. In light of these correspondences, we shall use the suggestive notation and terminology 
indicated in Table 2 for the various actions of the discrete diffeomorphism group on its Lie algebra, its Lie 
algebra's dual, and the space of discrete scalar fields. 
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Abstract formulation 


Continuous 


Discrete 


G 


Diff vol (M) 


V(M) 





X div (M) 


0(M) 


9* 


1 (M)/rfO°(M) 


f2j(M)/dn»(M) 


(■,■) :q* xg^K 


/ b / \ ? 

/ w b (v)<ix 
Jm 


Tr(C bT £lB) 


Ad g ?7 






Ad*-lyU 






ad^?7 


-£ u v 




ad^ 


£ u w b 


-[^.C^fi- 1 


G-action on scalar fields 







Table 1: Table of correspondences between the continuous and discrete diffeomorphism groups. 



Notation 


Meaning 


Terminology 


q*B 


qBq- 1 


Pushforward of a discrete vector field 


q*B 


q~ l Bq 


Pullback of a discrete vector field 


q*C b 


qC^flq^n- 1 


Pushforward of a discrete one-form 


q*C b 




Pullback of a discrete one-form 


£ A B 


-[A,B] 


Lie derivative of a discrete vector field 


£a& 




Lie derivative of a discrete one-form 


q*F 


qF 


Pushforward of a discrete zero-form 


q*F 


q-^F 


Pullback of a discrete zero-form 



Table 2: Table of notation and terminology for the actions of the discrete diffeomorphism group on its Lie 
algebra, its Lie algebra's dual, and the space of discrete scalar fields. 
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2.3 Nonholonomic Constraints 

Define the constrained set S C D(M) to be the set of matrices A E 0(M) whose entries are nonzero 
only if cells Ci and Cj are adjacent. In addition to enhanced sparsity, the set S has the physically appealing 
feature that its members encode exchanges of fluid particles between adjacent cells only. As Pavlov and 
co-authors [39] show, a nonzero entry Ay of a matrix A G S that discretizes a vector field v <E X<n v {M) 
represents (up to a constant multiplicative factor) the flux of the field v across the face shared by 
adjacent cells Ci and Cy. 



In applications, we will use matrices only in S to represent velocity fields, bearing in mind that this choice 
corresponds to the imposition of nonholonomic constraints on the dynamics of the system; indeed, the 
commutator of a pair of matrices in S is not necessarily an element of S. 

The imposition of nonholonomic constraints will be effected by constraining the variations in our deriva- 
tions of variational integrators; this, in turn, will amount to replacing equalities on discrete vector fields 
and one-forms with weak equalities in the sense laid forth below. These notions will all become clearer in 
Section 3 when we derive a family of integrators by taking variations of a discrete action and equating those 
variations to zero to obtain numerical update equations. These equations will later be replaced by weak 
equalities for the purposes of numerical implementation. 

2.3.1 The Discrete Flat Operator 

Prior to this section, we have employed the notation C b to refer to arbitrary elements of f2j(M)/df2^(M). We 
shall now introduce an explicit meaning for the symbol b by defining a discrete flat operator taking discrete 
vector fields to discrete one-forms. Designing a discrete flat operator in such a way that the discrete theory 
maintains its parallelism with the continuous theory is a nontrivial task and is perhaps one of the most 
important contributions of Pavlov and co-authors [39] toward the development of the discrete geometry of 
fluid flow. We give the definition of the discrete flat operator here, as well as an example of a flat operator on 
a cartesian mesh, but we refer the reader to [39] for its derivation and a generalization to irregular meshes. 

Definition 2.5. (Discrete Flat Operator.) Choose a parameter e that measures the resolution of a mesh 
M (e.g. the maximal length of a mesh edge). A discrete flat operator is a map b : S — > ilJ(M)/df^(M) 
taking discrete vector fields (satisfying the nonholonomic constraints ) to discrete one-forms that satisfies 



for any A, B,C G S that approximate continuous vector fields u, v,w, where the limits above are taken as 
the mesh resolution e tends to zero and the parametric dependence of A, B, C , and the flat operator on e have 
been denoted via subscripting. 

Example. On a regular two-dimensional regular cartesian grid with spacing e, the operator b : S — > 




(2.25) 



(C b e >,B e ) -> (w b ,v) 
{& t %£ A B t ) -> (w\£ u v 
(C\B) = (B\C) 



(2.26) 
(2.27) 
(2.28) 



f^(M)/df^(M) defined by 




2e 2 C l3 if j G N(i) 
]T (C ik +C kj ) i£j€N(N{i)) 



(2.29) 
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is a discrete flat operator, where w^j = 1 if cells Ci and Cj are share a single vertex and = 2 if cells 
and Cj belong to the same row or column. 

Remark. Throughout this paper, we continue to make no notational distinction between arbitrary elements 
of Q d (M) / dQ d (M) and flattened elements of S, even though the latter quantities constitute a proper subset 

of n\(M)/dn° d {M). 



2.3.2 Weak Equalities 

Throughout this report, our use of variational principles for the derivation of numerical integrators will 
lead to weak equalities on discrete forms and discrete vector fields. Such weak equalities will be denoted 
using the hat notation = to remind the reader of the absence of a strong equality. A discrete one-form 
C* b G n 1 d (M)/dn° d (M) will be called weakly null, i.e., 

& = 0, (2.30) 

iff is zero when paired with any vector field in the constraint space S, i.e., 

(C b ,Z)=0 VZeS. (2.31) 

The equality (2.31) holds not for all Z € t)(M) but for all Z in a subset of 0(M), the constraint space S, 
hence the weak equality. Similarly, a discrete vector field B £ D(M) will be called weakly null, i.e., 

B = 0, (2.32) 

iff B is orthogonal to every Z in the constraint space S, i.e., 

(Z b ,B)=0 yZeS. (2.33) 



Weak Equalities on Discrete One- Forms. The following lemma characterizes the nature of solutions 
to equalities of the form (2.30). 

Lemma 2.6. Suppose C b is a discrete one-form satisfying (C b , Z) = for every Z eS, i.e. C* b = 0. Then 
there exists a discrete zero-form P for which 

C$j =Pi~ Pj (2-34) 

for every pair of neighboring cells Ci , Cj . 

Proof. This follows directly from the quotient space structure of 0(M)* = fi^(M) / dCl d (M.) together with 
the sparsity structure of Z. M 



Weak Equalities on Discrete Vector Fields. In a similar vein, equalities of the form (2.32) require 
special care. Most relevant to our studies will be the situation in which the quantity B in (2.32) has the 
form B = X — A for a fixed matrix A 6 [S,S] and an undetermined matrix X € 2?(M). In such a scenario, 
the solution X E 2?(M) is in general not unique. While choosing X = A certainly ensures the satisfaction of 
the equality (2.32), we are often interested in solutions X that belong to S, the physically meaningful space 
of discrete vector fields satisfying the nonholonomic constraints. 
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Figure 2.1: Arrangement of cells with indices i, j, ki,k2,h,l2,h, Ia in definition (2.36) of the sparsity operator 
on a two-dimensional regular cartesian grid. The corresponding arrangement of cells for a pair of vertically 
adjacent cells Ci and Cj is given by the 90-degree counterclockwise rotation of the diagram above. Further 
rotations in increments of 90 degrees give the remaining two possible orientations of the pair of adjacent 
cells Ci and Cj . 

To achieve this goal, we define a sparsity operator to be a map ^ : [S, S] — >■ S satisfying 



for any A 6 [S,S], Z E S. Given such a sparsity operator, a physically meaningful solution to the matrix 
equation X — A = is given by X — A^ . Note that (up to a rescaling of nonzero matrix entries) the sparsity 
operator ^ is nothing more than the dual of the fiat operator with respect to the il- weighted Frobenius inner 
product. 

Example. On a regular two-dimensional regular cartesian grid with spacing e, the operator ^ : [S, S] — > S 
given by 



is a sparsity operator, where fci, fea, Zi, ?2> ^3) ^4 are the indices of the six cells of distance no more than two 
from both cell i and cell j, as depicted in Fig. 2.1. Notice that the sparsity operator merely accumulates a 
weighted sum of all two-away transfers that pierce the interface between cells Ci and Cj . 

2.4 Discrete Loops 

In deriving conservation laws (namely Kelvin's circulation theorem and its variants) for fluid flows, a special 
role is played by the double dual Xdiv(M)** of %div(M) and its discrete counterpart 0(M)**. In the continuous 



(Z b ,A) = {Z\A X ) 



(2.35) 




A klJ + A lk2 + \{A hj + A lh + A hj + A iU ) if j e N(i) 



otherwise 



(2.36) 
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case, we will frequently identify Xdi v (M)** with the space of closed loops in M via the pairing 

( 7 ,w b )=i'w b , (2.37) 

where 7 : S 1 -> M is a closed loop in M and w b G £ d iv(Af)* = Q 1 (M)/dQ°(M). 

In the discrete case, we simply identify the double dual of U(M) with itself. This is consistent with 
Arnold's [5] treatment of Kelvin's circulation theorem, which, roughly speaking, views the circulation of a 
one-form w b around a closed loop 7 as the limit of the pairing of w with a family of "narrow currents" 
{v e } C %div(M) of width e as e — » 0. Here, v e denotes a vector field tangent to 7 that vanishes outside a 
strip of width e containing 7. 

Throughout this paper, we use the term space of discrete loops when referring to U(M)** = 0(M). 

3 Euler Poincare Equations with Advected Parameters 

Having discussed the geometry of the diffeomorphism group Diff vo i (M) and its discrete counterpart, we now 
present the theory of Euler-Poincare reduction. This framework will later be used to precisely describe the 
manner in which solutions to the governing equations for ideal fluid flow may be realized as geodesies on 
Diffvoi (M), or, equivalently, as solutions to the Euler-Poincare equations on Diff vo i(M) with respect to a 
kinetic energy Lagrangian t—\ f M ||v|| 2 <ix. 

In anticipation of our eventual extension to MHD and complex fluid flow, we present a generalization of 
Euler-Poincare theory dealing with Euler-Poincare systems with advected parameters. We present this theory 
in its most general form, namely that in which one is dealing with a physical system whose configuration space 
is a Lie group G, whose advected variables belonging to the dual of a vector space V, and whose Lagrangian 
L : TG x V* — > K is left- or right-invariant. We will see in Section 4 that, with the appropriate choice of the 
group G, vector space V, and Lagrangian L, a variety of governing equations arising in continuum dynamics 
can be recast as Euler-Poincare equations with advected parameters, including the governing equations for 
MHD, nematic liquid crystal flow, and microstretch continua. 

In the latter half of this section, we derive a structure-preserving discretization of the Euler-Poincare 
equations with advection. Our derivation directly extends the work of Bou-Rabee and Marsden [12], who 
derived variational Lie group integrators for ordinary Euler-Poincare systems without advection. 

To maintain consistency with Bou-Rabee and Marsden [12], we shall perform derivations for systems with 
left-invariant Lagrangians whose configuration groups act from the left on the space of advected parameters. 
Note, however, that most examples of semidirect products arising in continuum mechanics (including MHD 
and complex fluid flow) involve systems with right-invariant Lagrangians whose configuration groups act 
from the right on the space of advected parameters. We shall take care to present both the left-left and the 
right-right versions of the theorems that follow, proving the left-left versions and leaving the proofs of the 
right-right versions as exercises for the reader. 

3.1 The Continuous Euler-Poincare Equations with Advected Parameters 

The following theorem is due to Holm and co-authors [27]: 

Theorem 3.1. Let G be a Lie group which acts from the left on a vector space V . Assume that the function 
L : TG x V* — > K is left G-invariant, so that upon fixing a *E V* , the corresponding function L ao : TG — > K 
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defined by L an (v g ) :— L(v g , do) is left G ao -invariant, where G ao C G denotes the isotropy group of a Q . Define 
£:qxV* -> K by 

Kg~ lv g, S -1 ao) = ao). (3.1) 
Let a(t) = g(t)~ l a$ and £(t) = g(t)~ l g{t). Then the following are equivalent: 
1. Hamilton's variational principle 

8 ( L ao (g(t),g(t))dt = (3.2) 



o 



holds for variations of g(t) vanishing at the endpoints. 

2. g(t) satisfies the Euler- Lagrange equations for L ao on G. 

3. The reduced action integral 



T 

e(£,(t),a(t))dt (3.3) 



o 

is stationary under variations of (£(t), a,(t)) of the form 

6Z = V+[Z,V\ (3-4) 
5a = —rja, (3-5) 

where rj{t) is an arbitrary curve in g vanishing at the endpoints. 

4- The Euler-Poincare equations hold on g x V* : 

d 51 5£ 51 

<«) 

where o : V x V* — > g* is the bilinear operator defined by 

(va, v) v , xV = -(voa, ?7) B *x B (3-8) 

for all v G V, a G V* , and r\ G 0. 
Remarks. 

1. If L is right invariant and G acts on V from the right, then equations (3.4-3.7) read 

S£ = V -[£,*] (3.9) 
5a = -a<q (3.10) 

% = (3.12) 

where = 3(i)<?(i) _1 , ait) = aog(i) -1 , and o : V x V"* — > g* is the bilinear operator defined by 

(a7?,u)v» x v = -{voa,ri} g * XSi (3.13) 

for all v G V, a G V*, and rj G g. 
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2. There are also corresponding versions for right-invariant Lagrangians with G acting on V from the left, 
and left-invariant Lagrangians with G acting on V from the right; for the purposes of this study, only 
the left-left and right-right cases will be of interest. 

3. As mentioned in Section 2.3, we will ultimately impose nonholonomic constraints on the dynamics of 
discrete vector fields and discrete one-forms in our studies of fluid flow. Such constraints can be effected 
via two equivalent means: (1) constraining the variations 77(f) to lie in the nonholonomic constraint 
space, or (2) replacing the equality (3.6) with a weak equality in the sense laid forth in Section 2.3.2. 
For a full development of nonholonomic variational principles, see the text of Bloch [7]. 

Proof. We only sketch the proof here; see Holm and co-authors [27] for a full exposition. 

The equivalence of (1) and (2) is well-known. To prove the equivalence of (1) and (3), one uses the 
fact that generic variations of g(t) induce non-generic variations at the level of the Lie algebra. For matrix 
groups, this is seen as follows: Letting r\ = g~ 1 (Sg) shows that fj = g~ 1 {b~g) — g~ gg^ x {5g) = g~ 1 {5g) — £»7 
and hence <5£ has the form 5£ — g~ 1 {5g) — g 1 {Sg)g^ 1 g = fj + — v£, — V + [€■> v}- 

Lastly, to prove the equivalence of (3) and (4), take variations of (3.3), substitute the relations (3.4-3.5), 
and integrate by parts. ■ 



3.2 The Continuous Kelvin-Noether Theorem 

The following theorem, again due to Holm and co-authors [27], describes a conservation law that may be 
viewed as an extension of Noether's theorem. 

Theorem 3.2. Let G be a group which acts from the left on a manifold C and suppose JC : C x V* — > g** is 
an equivariant map; that is, 

K(gc,ga)=M* g *K{c,a) (3.14) 

for any g € G, c G C, and a€ P, where the actions of G on C and V* have been denoted by concatenation. 
Suppose the curve (g(t),£(t),a(t)) satisfies the (left-left) Euler-P oincare equations (3.6-3.7) with an initial 
advected parameter value a £ V* . Fix c Q e C and define c(t) = g(t)~ l c . Define the (left-left) Kelvin- 
Noether quantity I 11 : C x g x V* — s- R by 

I ll (c,Z ) a) = (£(c,a) ) ^\. (3.15) 

Then the quantity I 11 (t) := I 11 (c(t) , , a(t)) satisfies 

d -l"{t) = (lC{c,a),^oa\. (3.16) 



dt \ 5a 

Similarly, suppose G acts from the right on a manifold C and the map K, : C x V* — > g is equivariant in the 
sense that 

fC(cg, ag) = Adg*i/C(c, a) (3-17) 

for all g G G, c G C, and a G V* . Suppose the curve (g(t),£(t),a(t)) satisfies the (right-right) Euler- 
Poincare equations (3.11-3.12) with an initial advected parameter value ao G V* . Fix cq G C and define 
c(t) = Cog(t)~ l . Define the (right-right) Kelvin-Noether quantity I rr : C x g x V* — > K by 



r r ( C ,C,a) = (/C( C ,a),-). (3.18) 
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Then the quantity I rr (t) := I (c(t) , , a(t)) satisfies 

> = (^^M < 3 " 19 > 

Proof. A proof of this theorem is given in the work of Holm and co-authors [27]; the proof uses the 
equivariance of IC together with the Euler-Poincare equations to differentiate /" and manipulate into 
the form of (3.16). ■ 

Corollary 3.3. In the absence of the advected quantity a(t), the Kelvin- Noether quantities I : C x g — > K 
and I rr :Cxg-)M defined analogously as above are conserved by the Euler-Poincare equations (3.6-3.7) 
and (3.11-3.12), respectively. 

Remarks. 

1. In the absence of advected parameters, the Kelvin-Noether theorem is an instance of Noether's theorem: 
it relates the symmetry of the Lagrangian under the left (respectively, right) action of the group G on 
itself to the conservation of a momentum Ad*-i || (respectively, Ad* ||). 

2. In the context of ideal fluid flow, we shall see that the Kelvin-Noether theorem gives rise to Kelvin's 
circulation theorem: the line integral of the velocity field along any closed loop moving passively with 
the fluid is constant in time. 



3.3 The Discrete Euler-Poincare Equations with Advected Parameters 

A discretization of the Euler-Poincare equations with advected parameters may be obtained via a small 
modification of the discrete reduction procedure detailed by Bou-Rabee [11]. Let t : g — > G be a local 
approximant to the exponential map on G (hereafter referred to as a group difference map), let h £ K be a 
time step, and define the reduced action sum s^° : G K+1 — > K with fixed oq £ V* by 

K-l 

s?({9k}fU)= E^-< a ^ ( 3 - 2 °) 

fc=0 

with 

&=t- 1 (j»Wi)/A (3-21) 

and 

a k = 9k la o- (3.22) 

Notice that (3.20) provides a direct approximation to the action integral (3.3) via numerical quadrature. 
Hamilton's principle in this discrete setting states that 

Ssj a = (3.23) 

for arbitrary variations of gu £ G subject to Sgo = Sgx = 0. 

Before proceeding with a derivation of the discrete update equations arising from this principle of sta- 
tionary action, let us recapitulate a few properties of group difference maps r : q — > G derived in [11]. For 
reference, the basic Lie group theoretic notations appearing below are fixed in Appendix A.l. 
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Definition 3.4. (Group Difference Map.) A local diffeomorphism t : g — > G taking a neighborhood Af of 
£ g to a neighborhood of the identity e £ G with r(0) = e and t(£) -1 = r(— £) /or all £ £ Af is called a 
group difference map. 

Definition 3.5. (Right-Trivialized Tangent.) The right-trivialized tangent dr : g x g —> g o/ a group 
difference map r is defined by the relation 

DT(£)-5 = TR T(0 dT S (6). (3.24) 

This definition decomposes the differential of r(£) into a map dr^ on the Lie algebra and a translation from 
the tangent space at the identity of G to the tangent space at r(£). The right-trivialized tangent satisfies the 
following useful identity: 

dr 6 {5) = M T[i) dT^{5). (3.25) 

Definition 3.6. (Inverse Right-Trivialized Tangent.) The inverse right-trivialized tangent dr^ 1 : g x 
g — > g of a group difference map t is defined by the relation 

Dr-\g) ■ S = dT-\TR T{ _ a 6) (3.26) 

with g — t(£). This definition decomposes the differential ofr~ 1 {g) into a translation from the tangent space 
at g to the tangent space at the identity of G and a map dr7 l on the Lie algebra. The inverse right-trivialized 
tangent satisfies the following useful identity: 

dT^(5) = dTl ? 1 (Ad r( _ 4)( 5). (3.27) 

A convenient formula for computing the inverse right-trivialized tangent of a given group difference map 
t is given by 



*T'ro = 1 



■'(expfeJMO). (3.28) 



e = 



For proofs of the properties enumerated above, the reader is referred to [11] and the references therein. 
We now present the main theorem of this section, which provides a discrete analog of the continuous Euler- 
Poincare principle with advection discussed in Section 3.1. 



K 



Theorem 3.7. Let G act from the left on the vector space V and assume that the function L : TG x V* 
is left-invariant. Let £ : g x V* — > K be the left-trivialization of L as in Theorem 3.1. Suppose {g k , £fc, Qfc}fc = o 
is a discrete path for which (3.20) is stationary under variations of g k £ G with 5 go — 5gx = and with 
£k and ak given by (3.21-3.22). Then the sequence {g k , a fc}fcLo satisfies the (left-left) discrete Euler- 
Poincare equations with an advected parameter: 

9k+i = 9kT(h^ k ) (3.29) 

afe+i = r(-h( k )ak (3.30) 

f\P AP AP 

(^lr Wk =(drZl ik J^ + h^a k . (3.31) 

Similarly, suppose G acts on V from the right, L is right-invariant, and and a k given by 

Z k = T- 1 (g k+1 gZ 1 )/h (3.32) 
O/t = aogZ 1 - (3.33) 
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Then the stationarity of (3.20) implies that the sequence {g k , a-k} k= o satisfies the (right-right) discrete 
Euler-Poincare equations with an advected parameter: 

9k+i = T(h£ k )g k (3.34) 
afc+i = a k r(-h£ k ) (3.35) 

(*-«.)'| < 3 - 36 > 

where £ is the right-trivialization of L. 

Notice that (3.30-3.31) and (3.35-3.36) are decoupled from the group update laws (3.29) and (3.34), 
respectively. 

Proof. We prove the left-left version of the theorem here; the corresponding proof for the right-right version 
is similar. 

Let us first compute the variations S^ k and 5a k induced by variations of g k , in accordance with rela- 
tions (3.21-3.22). In terms of the quantity n k :— g k 1 Sg k , 

Sa k = -g^Sgkg^ao 

= -VkO-k- (3-37) 

Similarly, 

= Di- 1 (T(h£ k )) ■ (-TR r{h(k) i lk + TL T ( hM r) k+ i)/h 

= L>t _1 (t(/i^)) ' (-TR T(Mk) i] k + TR T{h( . k )PL& T{hl . h) ri k+ i)/h. (3.38) 
Applying definition (3.26) together with the identity (3.27), the latter expression reduces to 

= -dr~l (Tft) + drZ l Hk (Vk+i ) ■ (3.39) 

We are now equipped to compute variations of the action (3.20) under variations of {g k } k=0 with fixed 
endpoints. In terms of the diamond operator defined by (3.8), 

K-l i <■„ i K-l 



I 5£ \ £ ^— ± / S£ \ 

Ss 7 = J2-{ % flfe ' h Ja( / + ^\5K , ~ dT ^ (%) + dTl >^ ) 
fc=0 ^ 1 fc=o ^ q 1 



fe=0 1 ' fe=0 

Using (discrete) integration by parts together with the boundary conditions t]q = t\k = gives 



(3.40) 



(3.41) 



This implies that 
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must hold for each k = 1,2, . . . , K — 1 in order to ensure the stationarity of the action s^°. This relation 
together with rearrangements of (3.21-3.22) constitute the discrete Euler-Poincare equations (3.29-3.31) with 
an advected parameter. ■ 

Remark. This theorem holds for nonholonomically constrained systems, modulo the following amendment: 
we restict the variations {f?fc}jfL to lie in the nonholonomic constraint space, leading to the replacement 
of (3.31) and (3.36) with weak equalities in the sense of Section 2.3.2. 



3.4 The Discrete Kelvin-Noether Theorem 

The following theorem describes a discrete analogue of the continuous Kelvin-Noether theorem (3.2). We 
make a slight departure from the letter of the continuous theory by identifying g** with g. This is done in 
anticipation of our eventual discretizations of continuum theories, where we perform temporal discretizations 
after performing spatial discretizations of the relevant configuration spaces; needless to say, the spatially 
discretized groups will be finite-dimensional, implying g** = g. 

Theorem 3.8. Let G act from the left on a manifold C and suppose JC : C x V* — » g is an equivariant map; 
that is, 

IC(gc, ga) = Ad 9 /C(c, a) (3.43) 

for any g G G, c G C, and a G V* , where the actions of G on C and V* have been denoted by concatenation. 
Fix a group difference map r : g — > G and suppose the sequence {gk, £&, o, k }^ =0 satisfies the (left-left) discrete 
Euler Poincare equations (3.29-3.31) with an initial advected parameter value oq G V* . Fix cq G C and 
define c k = g k 1 c a . Define the (left-left) Kelvin-Noether quantity I : C x g x V* ^ M. by 

I U (c,£,a) = ((dT^y^Mca))- (3-44) 

Then the quantity ij} := I u (ck,£,k,cik) satisfies 

Ik - 4-1 = (hj^ o a k , fC(c k , a k )^ . (3.45) 

Similarly, suppose G acts from the right on a manifold C and the map K, : C x V* — > g is equivariant in 
the sense that 

K.{cg, ag) = Ad ff -i /C(c, a) (3.46) 

for all g € G, c € C, and a G V* . Suppose the sequence {g k , a-k] k= Q satisfies the (right-right) discrete 
Euler Poincare equations (3.34-3.36) with an initial advected parameter value oq G V* . Fix cq G C and 
define c k — cog^ 1 . Define the (right-right) Kelvin-Noether quantity I : C x g x V* — > M. by 

I rr {c,i,a) = UdTZttT^M^Y (3-47) 

Then the quantity /£ r := I rr (c k , a k ) satisfies 

IV - IIU = U^<>a k ,K{c k ,a k )\ . (3.48) 
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Proof. Again, we prove the left-left version of the theorem here and remark that the right-right version is 
nearly identical. Since K, is equivariant, 

= (i dT hl k )*^^( c k> a k)j ~ ($^^' dT hz k -M d T(hik-iM c k,a k ))^ (3.49) 
by (3.29). Using the identity (3.27) allows us to write 

I l k l ~I l Li = ((^J*||,/C(c fe) a fc )\ - (^^ dT zl 6k _Snck,a k ))\ 

= ((^)*^,£(c fe ,a fc )) (fd T ^_J*^Mc k ,a k )\ . (3.50) 

The result (3.45) then follows from (3.31). ■ 

Corollary 3.9. In the absence of the adverted quantity a k , the Kelvin- Norther quantities I 11 : C x g — > K 
and I rr :Cxg^R defined analogously as above are conserved by the discrete update equations (3.29-3.31) 
and (3.34-3.36), respectively. 

Remark. A careful reading of the above proof shows that the validity of Theorem 3.8 persists for nonholo- 
nomically constrained systems whenever JC(c k , a k ) belongs to the prescribed nonholonomic constraint space. 
We shall invoke this observation tacitly throughout Section 4. 

3.5 Symplecticity of the Discrete Flow 

The following paragraphs prove, via a small extension of the argument presented in [11], the symplecticity 
of the flow of the discrete Euler-Poincare equations (3.29-3.31) with advected parameters. 

For fixed a € V* , let Ad denote the set of sequences {gk,£,k, a k} k= o satisfying the (left-left) discrete Euler 
Poincare equations (3.29-3.31). Let : Ad —> K denote the restriction of the reduced action sum (3.20) to 
sequences in Ad, and define the initial condition space to be the set 

I = {(9o, to) I 9o eG.foej}. (3.51) 

Since the update equations (3.29-3.31) define a discrete flow Fd : Z x 1 — > I, we may equivalently think of 
as a real-valued function : I — » E on the initial condition space. 

Referring back to the proof of Theorem 3.7, the differential of s a d ° then consists of the boundary terms 
that were dropped in passing from (3.40) to (3.41): 

ds a d ° ■ (<W£o) = /(drZ^ K _ i r i ^--,9] ( 1 S9 K \ - /(dr^T^g^Sgo') + (h-^oao,gQ l Sg \ . (3.52) 

Define the one-forms 0+, 9~ : Tl -> K by 

e+(g k ,Ck) ■ (Sg k ,S^ k ) = ((dT-lr^g-^g^ - (h-^oa k ,g^Sg k \ (3.53) 

e-(9k,Ck) ■ (Sg k ,8^ k ) = hdrZl 6 J*^,r(-hC k )g k \F^Sg k \ (3.54) 
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so that 

ds a d ° ■ (5g ,6£ ) = 0-{g K -i,ZK-i) ■ (6g K -i,5^ K -i) - + (g ^ Q ) • (Sg ,S^). (3.55) 
In terms of the discrete flow Fd, 

ds a d ° = {Ff^ye- - e + . (3.56) 

Differentiating this relation and using the fact that the pullback commutes with the exterior derivative, we 
find 

= d 2 s a d " = (Ff- x )*d0- - d0+. (3.57) 

Now notice that d9~ and d9 + are identical: Differentiating a single term of the discrete action sum and 
recalling the steps that led to (3.40) in the proof of Theorem 3.7, we have 

d(£(£ k ,a k )h) = -0-(g k Ak) + O + (g k ,Ck) (3.58) 

and hence d9~ = d9 + since d 2 = 0. Defining the two-form aid = d0~ = d6 + then shows that 

(Ff- 1 )*^ = u d . (3.59) 

In other words, the discrete flow preserves the discrete symplectic two-form u)d- A similar argument may be 
used to prove the symplecticity of the flow defined by the right-right discrete Euler-Poincare equations with 
advected parameters. 

4 Applications to Continuum Theories 

This section leverages the theory presented in Sections 2 and 3 in order to design structured discretizations of 
four continuum theories: ideal, incompressible fluid dynamics; ideal, incompressible magnetohydrodynamics; 
the dynamics of nematic liquid crystals; and the dynamics of microstretch continua . For each of these 
continuum theories, we present (1) the continuous equations of motion for the system, (2) the geometric, 
variational formulation of the governing equations, (3) an application of the Kelvin-Noether theorem to the 
system, (4) a spatial discretization of the governing equations, (5) a spatiotemporal discretization of the 
governing equations, and (6) an application of the discrete Kelvin-Noether theorem to the discrete-time, 
discrete-space system. 

The continuous variational formulations that we present here for ideal fluid flow and MHD are well- 
documented; see references [5, 26, 27, 35] for a comprehensive exposition. The variational formulation of 
nematic liquid crystal flow and microstretch fluid flow is a more recent development due to Gay-Balmaz & 
Ratiu [21]. 

The reader will discover, upon completion of a first reading of this section, that the procedure employed 
for each of our discretizations adheres to a canonical prescription. While the details of this prescription will 
be made clearer through examples, we summarize here for later reference our systematic approach to the 
design of structured discretizations of continuum theories: 

1. Identify the configuration group and construct a finite-dimensional approximation G to that group. 
In the applications we present, this typically involves replacing the diffeomorphism group Diff vo i(Af) 
with the matrix group 2?(M) and replacing scalar fields with discrete zero-forms. 

x As mentioned earlier, our numerical treatment of the last two examples of continuum theories will only be valid in 2D; 
extension to three and higher dimensions requires an alteration of our discretization of Lie-algebra-valued fields that will be 
addressed in a different paper. 
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2. Compute the infinitesimal adjoint and coadjoint actions on g and g* , respectively. 

3. Identify the space of advected parameters and construct a finite-dimensional approximation V* to that 
space. Design an appropriate representation of G on V*. 

4. Compute the diamond operation associated with the action of GonV*. 

5. Write down the discrete-space Lagrangian I : g x V* — s- K. 

6. Compute the discrete-space, continuous-time Euler-Poincare equations (3.11-3.12) with advected pa- 



If a variational temporal discretization is desired, one continues in the following manner: 

7. Compute the exponential map exp : g — > G. If G = 2?(M), this is simply the usual matrix exponential. 

8. Construct a local approximant r : q — > G to the exponential map. It is often convenient to invoke the 
Cayley transform for this purpose. 

9. Compute the inverse right-trivialized tangent dr _1 : g x g —> R of t, as well as its dual. 

10. Write down the discrete-space, discrete-time Euler-Poincare equations (3.34-3.36) with advected pa- 



For finite-dimensional systems, only the last four steps are relevant. To illustrate the process of temporal 
discretization described in the last four steps, we will begin this section by presenting two finite-dimensional 
mechanical examples, the rigid body and the heavy top. These mechanical examples are not merely peda- 
gogical. Quite remarkably, the discrete equations of motion for the rigid body will be seen to exhibit strong 
similarities to the discrete fluid equations expressed in our framework. Another parallel between the heavy 
top and MHD will be highlighted. These analogies are manifestations of the fact that the governing equa- 
tions for perfect fluid flow and those for the rigid body are both Euler-Poincare equations associated with a 
quadratic Lagrangian, albeit on a different Lie group. The transparency of this correspondence in the discrete 
setting is one of the unifying themes of this section, and it serves as a strong testament to the scope and 
power of geometric mechanics in numerical applications. 

4.1 Finite-Dimensional Examples 
4.1.1 The Rigid Body 

Consider a rigid body moving freely in empty space. In a center-of-mass frame, its configuration (relative 
to a reference configuration) is described by a rotation matrix R G SO (3). Its angular velocity in the body 
frame is related to its configuration according to the relation 



where ft — (l^, f2 2 , ^3) G K 3 and * : (R 3 , x) — > (so(3), [•, •]) is the Lie algebra isomorphism associating 
vectors in K 3 to antisymmetric matrices: 



rameters. 



rameters. 



Cl = R^R, 



(4.1) 







^3 




n = 



n 3 



o 



(4.2) 



fix 



/ 
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The rigid body Lagrangian is 

£(fl) = -Ml • n = ^Tr(f2 T jri), (4.3) 

where I 6 Sym(3) + is the inertia tensor of the body in a principal axis frame, and J is related to I via 

J = Tr(I)7 - I. (4.4) 

This Lagrangian is left-invariant since the map R h- > QR leaves Jl = R _1 R invariant for any Q € 50(3). 
Identifying so(3) with its dual via the pairing (v, w) = v • w = |Tr(v T w), define 

St - 

n = = nj + jn = in. (4.5) 

The Euler-Poincare equations (3.29) (without advected parameters) on so (3)* read 

^-[tl,Cl}=0, (4.6) 

or, equivalently, 

™ - n x n = 0. (4.7) 

Temporal Discretization. For a given group difference map r : so(3) — > 50(3), the discrete-time Euler- 
Poincare equations (3.29-3.31) (again without advected parameters) read 

R fc+ i = K k T(hn k ) (4.8) 
(Klrtl k = (dr : K ki ril k ^. (4.9) 

One choice for a group difference map t is the Cayley transform 

whose inverse right-trivialized tangent takes the form (see [11]) 

dcay^(#) = ( I - "J # f 1+ = * - *] - (4.11) 
Equation (4.9) with r = cay can then be written 

rV-ri^ _ [ft fc _ 1 ,n fc - 1 K[ft fc ,n fc ] + _ j = . (4 .i 2) 
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The Discrete Kelvin-Noether Theorem. To apply theorem (3.8), let G — SO(3) act on C := so(3) 
via the adjoint representation and let JC : so(3) — » so(3) be the equivariant map W h-> W. The discrete 
Kelvin-Noether theorem (free of advected parameters) then implies that 

<( dT ^)*n*,W fc ) (4.13) 

is independent of k for any Wo € so(3), where Wfc = R7 WoRfe- Rearranging (4.13) results in 

(Rk(dT-£ )*n fc Rfc\W >. (4.14) 



Since W is arbitrary, this shows that the rigid body update equations preserve the quantity 

jt k = R k (drJ- h yu k R-j- 1 , (4.15) 

which is a discrete analogue of spatial angular momentum of the body. In other words, the rigid body 
discrete update equations (4.8-4.9) preserve a discrete spatial angular momentum exactly. 

4.1.2 The Heavy Top 

The heavy top is a mechanical system consisting of a rigid body with a fixed point rotating in the presence 
of a uniform gravitational field in the — e3 direction. Let l\ denote the vector pointing from the top's 
point of fixture to its center of mass, where x is a un it vector and I is a scalar. Let R € SO (3) denote 
the configuration of the top relative to the upright configuration and let Jl denote the top's body angular 
velocity, related to R via Cl — R~ 1 R. 

If r € K 3 is the direction of gravity as viewed in the body frame, then T is a purely advected quantity: 

r(t)=K(t)- 1 r = R(t)- 1 e 3 . (4.16) 

We therefore treat T as an element of a representation space V* := M 3 , on which G := 50(3) acts from the 
left by multiplication: 

r r = Rr. (4.17) 

The infinitesimal action of so(3) on V* is then 

£i r = n x r. (4.18) 

Now note that for any wet 3 , 

(nxr)-w = -(wxr)-n, (4.19) 

showing that the diamond operation o : M 3 x M 3 —> so (3) (see 3.8) is given by 

woT = (w~xT) = [w,f]. (4.20) 

The Lagrangian £ : so (3) xM 3 4K for the heavy top is the kinetic energy of the top minus its potential 
energy: 

*(n,r) = i<in,n)-M 5 ir.x, (4.21) 

where M the mass of the top, g is the acceleration due to gravity, and I is the top's inertia tensor. 
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The Eulcr-Poincare equations (3.29-3.30) with the advected parameter T now read 

dii 



dt 



[U,n} = Mgl[T,x] 



dT 

~dJ 



or, equivalently, 



dU 

~dt 



- n x ft = MglT x x 



dT 

~d~t 



rxn. 



(4.22) 
(4.23) 

(4.24) 
(4.25) 



Temporal Discretization. The discrete-time Euler-Poincare equations (3.29-3.31) for the heavy top read 

R fc+1 = R fc r(/iO fe ) (4.26) 



r fc+ i = T(-hn k )T k 
-i \*tt. _ r^-i \* n 



(4.27) 
(4.28) 



hMgl[T k , X ]- 

To facilitate an eventual analogy with MHD, let us recast the second equation in terms of T € so(3): 

f fc+i = T(-htl k )t k T(hn k ) (4.29) 
Now use the Cayley group difference map (4.10) and rearrange (4.28-4.29) to obtain 

rife-rV! [n fc _ ll nj k _i] + [n fc ,n fc ] , h. 



Tfe — Tfe_i 
h 



Tfe-i + T k 



h 



-(njfe-injk-infc-i - n k u k n k ) = M g i[r k , x } (4.30) 



(n^Tk^n^ - n k ^f k n k ^) = o. (4.31) 



The Discrete Kelvin-Noether Theorem. To apply theorem (3.8), again let G — 50(3) act on C := 



so (3) via the adjoint representation and let K : so (3) 
The discrete Kelvin-Noether theorem then says that 



so(3) be the equivariant map (W, T) H> W. 



((dT-} k )*Tl k ,W k ) - ((dr-^yiLk-uWk-!) = hMgl{[t k , X },W k ), (4.32) 

where Wfc = R^ WoRfc = R7T Wo. That is, W k € M 3 is the body frame representation of a fixed 
spatial vector W G M 3 . Choosing W = T so that Wj. = T k causes the right-hand side of (4.32) to 
vanish, showing that ((dT~^ )*Tl k ,T k ) = (Rfc^T"^ ^n^R^ 1 , e 3 ) is constant. Thus, the discrete update 



hCl 



equations (4-26-4-28) preserve the e%- component of the heavy top's spatial angular momentum 



7r* = R*(dr-i )*n fc Ri 



(4.33) 
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4.2 Discretization of Continuum Theories 

We now proceed with our discretization of four continuum theories: ideal, incompressible fluid dynamics; 
ideal, incompressible magnetohydrodynamics; the dynamics of nematic liquid crystals; and the dynamics of 
microstretch continua. 

4.2.1 Ideal, Incompressible Fluid Flow 

The Continuous Equations of Motion. The governing equations for the dynamics of an ideal, incom- 
pressible fluid occupying a domain M C R 3 are the Euler equations: 

~ + (v • V)v = -Vp (4.34) 
V • v = (4.35) 

Equation (4.34) describes the temporal evolution of the fluid velocity field v in terms of the coevolving scalar 
pressure field p. The latter quantity is uniquely determined (up to a constant) from the constraint (4.35), 
which expresses mathematically the incompressibility of the fluid. 

Geometric Formulation. The configuration space of the ideal, incompressible fluid is the group G = 
Diff vo i(M). In terms of material particle labels X and fixed spatial locations x <E M, a passive fluid particle 
traces a path x(X, t) — ^t(X) under a motion ip t in G. The spatial velocity field v(x, t) is related to ipt via 

v(x,t)=ip f (i P t 1 (x)). (4.36) 

Recognizing (4.36) as the right translation of a tangent vector ip t € T Vt G to the identity, we regard v(x,t) 
as an element of g = X ( a v (M), the Lie algebra of G. 

The fluid Lagrangian I : g — > R, regarded as the right-trivialization of a Lagrangian L : TG — > R, is the 
fluid's total kinetic energy 

^(v) = i(v b ,v), (4.37) 

where the pairing above is that given by (2.4). 

The Euler-Poincare equations (3.11) (without advected parameters) read 

+ £ V ^ - -dp. (4.38) 

The pressure differential appearing here is an explicit reminder that the Euler-Poincare equations formally 
govern motion in g* = Q 1 (M)/dQ°(M), a quotient space in this case. In the language of vector calcu- 
lus, (4.38) has the form of equation (4.34), with p related to p by p = p + • v. 

The Continuous Kelvin-Noether Theorem. The ideal, incompressible fluid provides a canonical il- 
lustration of the utility of the Kelvin-Noether Theorem (3.2) in the absence of advected parameters. Take 
C to be the space of loops (closed curves) in the fluid domain M, acted upon by G from the right via the 
pullback, and take K, : C — » g** to be the equivariant map defined by 

(/C( 7 ),w b > = I w\ (4.39) 
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where 7 : S 1 — > M is a loop and w is a one-form. The Kelvin-Noether theorem then says that 

v b (t)=0 (4.40) 

(*) 

where v(t) is the fluid velocity field and j(t) is a loop advected with the flow. This is known classically as 
Kelvin's circulation theorem: the line integral of the velocity field along any closed loop moving passively 
with the fluid is constant in time. 




Spatial Discretization. To discretize (4.38), cover M with an mesh M and replace Diff vo i(M) with the 
finite-dimensional matrix group G = T>(M). Denote the discrete counterpart of ip by y and that of v by Y; the 
quantities y and Y are elements of the finite-dimensional diffeomorphism group and its Lie algebra q = 0(M), 
respectively. As such, y is an fi-orthogonal, signed stochastic matrix, and Y is an f2-antisymmetric, row-null 
matrix, where f2 is the diagonal matrix of cell volumes described in Section 2. 

The spatially discretized fluid Lagrangian I : g — > R is the fluid's total kinetic energy 

i{Y) = \{Y\Y), (4.41) 

where the pairing above is that given by (2.16). This Lagrangian is right invariant: in analogy with (4.36), 
the relation between y and Y is given explicitly by 

Y = yy-\ (4.42) 

showing that the map y 1— ?► yq with q 6 2?(M) leaves (4.41) invariant. 

The discrete-space, continuous-time Euler-Poincare equations (3.11) (without advected parameters) read 

— + £ Y Y^0, (4.43) 
ot 

where £yF k is the f2-weighted matrix commutator (2.20) and the hatted equality denotes a weak equality 
in the sense of Section 2.3.2. Namely, there exists a discrete scalar field P for which 

(^ + £ Y Y^ =-(dP) y (4.44) 

for every pair of neighboring cells Ci and Cj . 

Note the similarity between (4.43) and the rigid body equation (4.6) upon writing £yY b in its explicit 

form [y^y]^- 1 . 



Temporal Discretization. The discrete-space, discrete-time Euler-Poincare equations (3.34-3.36) (again 
without advected parameters) read 

Vk+i = T(hY k )y k (4.45) 
[drZlyjYt = (dr^TY^. (4.46) 

A convenient choice for a group difference map r : c)(M) — > 2?(M) is once again the Cay ley transform 

r(y) = cay(r) =(/-■£) '(^l)- ( 4 - 47 ) 
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It is well-known that the Cayley transform is a local approximant to the matrix exponential and preserves 
structure for quadratic Lie groups; that is, if Y is f2-antisymmetric, then cay(y) is f2-orthogonal. Conve- 
niently, the Cayley transform maps row-null matrices to signed stochastic matrices, making it a genuine map 
from the Lie algebra 0(M) to the Lie group 2?(M). To see this, note that if Y is row-null, then (i — y) 
and (/ + y) are each stochastic; now use the fact that the set of signed stochastic matrices is closed under 
multiplication and inversion. 

As we saw with the rigid body, the inverse right-trivialized tangent of the Cayley transform takes the 
form (see [11]) 

Y\ „ / Y\ / 1 „ \ „ 1. 



dcsyy^Z) = II - -j Z U + - 1 = 17 + -£ Y j Z - -YZY (4.48) 
The dual of efcay -1 with respect to the pairing (2.16) can be computed as follows: For any Z £ 0(M), X £ 

(X\dc aj y\Z)) = -Tr (tl (j +j) X 1 

-Tr I fiz 1 1 + — ] x b n ( i - — ) rr 1 



This shows that 



2 J V 2 

4i- Y - 

2 \ 2 



i+~)x*n(i-~)n- l ,z). (4.49) 



(dcay^ 1 )*^) = ( 1+ X b Q (i - CT 1 = (i - ^£ Y ^j X b - ^YX^QYQ- 1 . (4.50) 

Substituting into (4.46) and rearranging leads to the following discrete update equation for ideal, incom- 
pressible fluid flow: 

V b — Vl> /V Y* 4- /V Y} h 

k l k - 1 + Yk -> k -\ Yk k + fa^YUar^n- 1 Yt.Y^n- 1 ) = o. (4.5i) 

Notice once again the evident parallel between (4.51) and the rigid body update equation (4.12). 

The Discrete Kelvin-Noether Theorem. Let C be the space of discrete loops in M and let G = X>(M) 
act on C = g via the discrete pullback: T ■ y := y*T = y _1 Ty for any T £ C, y £ G. Let JC : C — > g be the 
equivariant map T t— > T. The discrete Kelvin-Noether theorem then says that the quantity 

M fc = ((drr£ n )*n b ,r*) (4.52) 

satisfies 

M k - M fe _! = (4.53) 

with Ft = r -y^ 1 = (yfc)*r , i.e. F^ is a discrete loop advected passively by the fluid flow. Identifying (4.52) 
as the circulation of the one- form (drZhy k )*~Yk along the curve r&, we see that (4.53) gives the discrete 
analogue of Kelvin's circulation theorem. In other words, a discrete Kelvin's circulation theorem holds 
for the discrete-space, discrete-time fluid equations laid forth by Pavlov and co-authors [39]. Notice the 
beautiful parallel between the discrete Kelvin's circulation theorem for the ideal fluid and the discrete 
angular momentum conservation law derived earlier (in Section 4.1.1) for the rigid body. 
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4.2.2 Ideal, Incompressible Magnetohydrodynamics 

The Continuous Equations of Motion. The governing equations for ideal, incompressible MHD, which 
studies the motion of a perfectly conducting incompressible fluid occupying a domain M C K 3 in the presence 
of a coevolving magnetic field B, are the ideal MHD equations [23]: 



dB 



dv 

— + (v • V)v = (V x B) x B - Vp (4.54) 
at 

- V x (v x B) = (4.55) 

V • B = (4.56) 

V • v = (4.57) 

Equation (4.54) is the classic Euler equation for the fluid velocity field v with an additional term (V x B) x B 
corresponding to the Lorentz force j x B acting on mobile charges in the fluid. The second equation (4.55) 
describes the evolution of the magnetic field B and has a particularly simple interpretation when recast in 
the form 

f)B 

— + £ V B = 0; (4.58) 

that is, the magnetic field is advected with the fluid flow. The constraints (4.56-4.57) represent the absence 
of magnetic charge and the incompressibility of the fluid flow. The latter constraint uniquely determines 
the pressure p appearing in (4.54), while the former constraint automatically holds for all time if it holds 
initially, as can be seen by taking the divergence of equation (4.55) and using the fact that V • (V x u) = 
for any vector field u. 

One may check through differentiation that the quantities 



M 



-V V + -B B ) <lx. (4.:.!)) 



and 



J = / v-Bdx (4.60) 

J M 

are constants of motion for the ideal MHD equations (4.54-4.57). The former quantity is the total energy 
of the system, being the volume-integrated sum of the kinetic energy density • v of the fluid and the 
potential energy density |B ■ B of the magnetic field. The quantity J is known as the cross-helicity of the 
fluid, which may be shown to bear a relation to the topological linking of the magnetic field B and the fluid 
vorticity V x v [5]. 



Geometric Formulation. As with the ideal, incompressible fluid, the configuration space for the ideal, 
incompressible magnetofluid is the group G = Diff vo i(M). Curves <pt in G encode fluid motions as usual, 
and the fluid's spatial velocity field v belongs to q — X^ V {M). 

According to (4.58), the magnetic field B is an advected parameter. With this in mind, we choose the 
representation space V = fl, (M)/d£l°(M) and identify V* with Xdiv(M) under the usual pairing. Elements 
(p of G act on V from the right via the pullback, so that the induced right action of G on V* is again given 
by the pullback: 

B-ip = tp*B. (4.61) 
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The infinitesimal action of g on V* is via the Lie derivative: 

B v = „£ v B. (4.62) 

The diamond operation (3.13) may be computed as follows: For any C b £ V, B g V* , and v£g. 

(B • v, C b )v*xv = (C b , £ v B) B * xg 
= -(C b , 4v) r x B 
= (,£ B C b , v) B . X0 

= -(C b oB,v) rxfl . (4.63) 

Hence, 

C b oB = -4C l . (4.64) 

The fluid Lagrangian I : g x V* — > K is the fluid's total kinetic energy minus the potential energy stored 
in the magnetic field: 

*(v,B) = I(v b ,v)-i(B b ,B), (4.65) 

where the pairing above is that given by (2.4). 

The Euler-Poincare equations (3.11) with the advected parameter B read 

dv b 

— + £ v v b = £ B B b - dp (4.66) 

— + £ V B = 0. (4.67) 

In the language of vector calculus, (4.66-4.67) have the form of equations (4.54-4.55), with p related to p by 
p = p+|v-v-B-B. 

The Continuous Kelvin-Noether Theorem. As with the ideal fluid, take C to be the space of loops in 
M and let G = Diff vo i(M) act on C via the pullback. Let JC : C x V* — s- q** be the equivariant map defined 

by 

(/C( 7 ,B),w b > = /w b , (4.68) 
Ji 

where 7 : S 1 — >■ M is a loop, B € V* , and w b is a one-form. The Kelvin-Noether theorem then gives the 
statement 

^ / v b (t) - / £ B B b (t) (4.69) 
at J l{t) J l{t) 

where v(t) is the fluid velocity field, B(t) is the magnetic field, and j(t) is a loop advected with the flow. 

A second momentum evolution law is obtained from the Kelvin-Noether theorem if one identifies q** with 
g = Xdiv(M) and takes C = £ < n v (M) and /C(u,B) = u, with G acting on Xdi v (-^0 via the pullback. The 
Kelvin-Noether theorem then gives 

|(v b (i),u(t)) = (£ B B b (<),u(0), (4.70) 
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where u(t) = ip(t) *u(0) is a vector field advected with the flow. Choosing u(0) = B(0) leads to the conclusion 

|(v b (i),B(t))H^ B B^),B(i)) 

= -(B b (i),£ B B(t)) 

= 0; (4.71) 

that is, the quantity 

J(t) = (v\t),B(t)) (4-72) 

is conserved along solutions to the governing equations for ideal MHD. As mentioned in this section's intro- 
duction, the quantity J is known as the cross-helicity in MHD, and can be shown to bear a relation to the 
topological linking of the vorticity field V x v and the magnetic field B [5] . 

Spatial Discretization. A spatial discretization of the ideal, incompressible MHD equations is obtained 
by replacing Diff vo i(Af) with the finite dimensional matrix group G — T>(M). As before, denote the discrete 
counterpart of ip by y and that of v by Y . Replace the representation space Jl 1 (M) / 'dil° '(M) and its dual 
X div (Af) with V = Q^(M)/dQ^(M) and V* = 0(M), respectively. Denote the discrete counterpart of B by 
R; the quantity R is an ^-antisymmetric, row-null matrix belonging to V* . 

Elements y of G act on V from the right via the discrete pullback, so that the induced right action of G 
on V* is again given by the discrete pullback: 

R-y = y*R = y~ 1 Ry. (4.73) 

The infinitesimal action of q on V* is via the discrete Lie derivative: 

R-Y = £ Y R=-[Y,R]. (4.74) 

The derivation of the continuous diamond operation (4.64) carries over seamlessly to the discrete setting, 
giving 

S b oi? = -£ R S\ (4.75) 

The spatially discretized MHD Lagrangian £ : g x V* —> R is the fluid's total kinetic energy minus the 
potential energy stored in the magnetic field: 

t(Y,R)= \{Y\Y)-\{R\R), (4.76) 

where the pairing above is that given by (2.16). 

The discrete-space, continuous-time Euler-Poincare equations (3.11) with the advected parameter R read 

— + £ Y Y b = £ R R b (4.77) 
C ^ + £ Y R = 0, (4.78) 



where P is discrete scalar field. 
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Temporal Discretization. The discrete-space, discrete-time Euler-Poincare equations (3.34-3.36) with 
the advected parameter R read 

Z/ fe+ i = T(hY k )y k (4.79) 
R k+1 =T(hY k )R k r{-hY k ) (4.80) 
{drZ l hYk )*Yl = {dT^Y kl )*Yl_ l + h£ Rk R\. (4.81) 

Using the Cayley transform (4.47) as a group difference map, the update equations (4.80-4.81) become, 
upon rearrangement, 

Y * + £y ^ Y ^ + £y ^ + ^ (Y k _ i Y^_ 1 f2F s ._ 1 J7 _1 - Y k Y k b nY k n-i) = £ Rk Rl (4.82) 

Rk '^ k - 1 +£y k ^ 1 f^l±3t\ + ^(Y k _ 1 R k _ 1 Y k _ 1 -Y k ^R k Y k _ 1 )=Q. (4.83) 

Referring back to equations (4.30-4.31), the reader can note the striking similarity between the heavy top 
update equations and the MHD update equations above. 

The Discrete Kelvin-Noether Theorem. Take C to be the space of discrete loops in M and let G = 

V(M) act on C = 5(M) via the pullback. Let K, : C x V* q be the equivariant map (T.A) ^ T. The 
discrete Kelvin-Noether theorem then says that the quantity 

Jk = ((dTll Yh yY k \T k ) (4.84) 

satisfies 

J k -J k - 1 = h{£ Rh R k ,T k ) (4.85) 
with T fc = T • y^ 1 = {y k )*T . If r = i? , then T k = R k and we find 

Jk — Jk-i = h(£ Rk R k ,R k ) 
= -h(R b k , £ Rk R k ) 

= 0. (4.86) 

Notice that (4.84) serves as the discrete analogue of the cross- helicity J = (v b ,B) arising in MHD. Thus, we 
have proven that the cross-helicity (4-84) is preserved exactly along solutions to the discrete-space, discrete- 
time MHD update equations (4-80-4-81). 



4.2.3 Nematic Liquid Crystals 

The dynamics of complex fluids differ from those of ordinary fluids through their dependence upon micro- 
motions, i.e., internal changes in the shape and orientation of individual fluid particles that couple with their 
ordinary translational motion. In the particular case of nematic liquid crystals, the fluid particles are rod-like 
and the internal motions are purely rotational. In a three-dimensional domain M 1 C IR 3 the variables of 
interest are the fluid velocity field v G 3L<i\ v {M'), the local angular velocity field v € T(M' ,1R 3 ), and the 
director n € F(M' , ]R 3 ). The latter variable describes the local orientation of fluid particles and its value at 
any point x € M' is restricted to have unit length, so that n may equivalently be viewed as a map from the 
domain to S 2 . 
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The Continuous Equations of Motion. In terms of these variables, the equations of motion for three- 
dimensional, incompressible, homogeneous nematic liquid crystal flow are given by 

| + (v.V)v = ^(^.Vn)-Vp (4.87) 

+ (v • V> = h x n (4.88) 

(Vn)v + n x v = (4.89) 
V • v = 0, (4.90) 



dn 



where p is a scalar pressure field, 



and F(n, Vn) is a scalar function describing the potential energy stored by irregularities in the fluid particles' 
alignment, hereafter referred to as the free energy. The notation Vn denotes the 3x3 matrix whose i th row 
is Vn^. 

The governing equations for nematic liquid crystal flow have the following interpretation: Equation (4.87) 
is the classic Euler equation for the fluid velocity v with an additional term accounting for irregularities in 
the director field that react back on the flow; equations (4.88-4.89) describe the advection of the local angular 
velocity field v and director field n; and (4.90) expresses the incompressibility of the flow. 

We will simplify the analysis considerably by restricting the rest of our discussion to two-dimensional 
flows with the director everywhere tangent to the plane. For concreteness, choose M C M' to be the x-y 
plane and define the variables u> <E F(M), a € F(M, S 1 ) via 

v = (0,0, cuf (4.92) 
n = (cosa,sina,0) T , (4.93) 

where we identify 5" 1 with the real numbers modulo 2ir. 

In terms of these variables, the equations of motion for two-dimensional, incompressible, homogeneous 
nematic liquid crystal flow are given by 

<9v / dF \ 

+ (v • V)v = -ft — Va - Vp (4.94) 



dt \da t i/ 
da 



dt 



v-Va = w (4.96) 
V • v = 0, (4.97) 



Remark. We will deal only with the two-dimensional case for the remainder of this section. Note that for 
higher dimensions, Lie algebra valued functions such as v will require a different discretization than the one 
used here in order to respect the group product involved in this model; we will explore this extension in a 
future paper. 
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Geometric Formulation. Rather remarkably, the governing equations (4.87-4.90) for the dynamics of 
nematic liquid crystals arise from a variational principle on a certain Lie group with respect to a standard 
Lagrangian given by the difference between kinetic and potential energies. Gay-Balmaz & Ratiu [21] show, in 
particular, that the appropriate configuration group for three-dimensional nematic liquid cyrstal flow is a con- 
tinuum analogue of the special Euclidean group, namely the semidirect product Diff vo i(M') ® F(M ' , SO (3)), 
where F(M' , 50(3)) is the space of maps from the fluid domain M' C M 3 to 50(3). 

The configuration space for two-dimensional nematic liquid crystal flow is the semidirect product 

G = ~DiK vol (M)®F(M,S 1 ), (4.98) 

where ^(Af^ 1 ) is the space of maps from the two-dimensional domain to 5 1 . The group product is given 

by 

fai,0i)(¥>a,02) = faiop 2 ,^0i+02) far(¥>i,0i) s (V2,02)eG. (4.99) 
The Lie algebra of G is 

g = X div (M)®T(M). 

We identify the space dual to g with 

g* = n 1 {M)/dn Q {M) x F(M) 

through the pairing 

((w b ,7T),(v>)) = (w b ,v) + M) 

for w b e n 1 (M)/dQ°(M), v e Xdiv(M), and tt, V G ?{M), where 

(w b ,v) = f w b (v)dx 

and 

(tt, ip) = / 7r(x)^)(x) dx. 

The bracket on g is computed to be 

ad( U;W ) (v, ip) = (-£ u v, £ v oj - £ u ijj), (4.105) 

with dual 

ad( UA)) (w b ,7r) = (£ u w b + vrdo;, £ u tt). (4.106) 

The space of advected parameters (not to be confused with the second factor of G) is T{M, 5 1 ) , whose 
elements a are acted upon from the right by elements (ip, 9) of G in the following manner: 

a-{ip,9)=cp*a-9. (4.107) 

One checks that (4.107) is a well-defined right action that is consistent with the group multiplication 
law (4.99). 

Note that the theory presented in Section 3 does not, strictly speaking, apply to this situation, as 
!F(M, 5 1 ) is not a vector space. Rather than introducing a flurry of new notation and terminology, we will 
simply regard members of T(M, 5 1 ) as real- valued functions on M when necessary. Those readers interested 



(4.100) 

(4.101) 
(4.102) 

(4.103) 
(4.104) 
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in the theoretical framework underpinning Eulcr-Poincare systems with nonlinear advected parameter spaces 
are encouraged to consult the work of Gay-Balmaz & Tronci [22]. 
The infinitesimal action of g on T(M, S 1 ) is given by 

a- (u,w) = £ u a-uj for (u, u) G X div (M) ® F(M). (4.108) 

The diamond operation (3.13) is computed as follows: For any a, j3 £ F(M 1 S 1 ) and (u,o>) € g, 

(a-(u,w),0)=(£ n a,0)-(w,0) 
= (pda,u) - 

= -((-(3da,(3),(u,oj)) (4.109) 

implying that 

/3oa= (-/3da,(3). (4.110) 

The reader may find it helpful at this point to refer to [27] for more examples of calculations of this type. 
The Lagrangian I : [Xa v (M) (S) F(M)] x F{M, S 1 ) — > R for nematic liquid crystals is 



11/" 

£(v,uj,a) = -(v b ,v) + ~(u],w) - / F(a,da)dx, (4.111) 



2 N ' ' 2 



where -F(a, da) is the free energy. The first two terms in the Lagrangian may be recognized as the kinetic 
energy of the fluid due to translational motion and internal rotations, respectively. 
The Euler-Poincare equations (3.11) with the advected parameter a read 

^— + £ v v b = ~?-da-dp (4.112) 
at da 

% + £*> = £ (4.113) 
at da 

da 

— + £ v a = u, (4.114) 
at 

where the term uidui — d(u> 2 /2) in (4.112) has been absorbed into the pressure differential dp. 

If (ip, 9) € G is used to denote the fluid configuration, then the evolution of a is given explicitly by 

a = a ■ (if., 9)- 1 = a Q ■ {ip' 1 , -tpj) = ^(a + 6). (4.115) 

In particular, if ao = 0, a and 9 are related via a transformation between spatial and material coordinates: 

aoip = 9. (4.116) 

A common choice for the free energy function F is 

F(a,da) = ^||cM| 2 , (4.117) 



which corresponds to the so- named "one-constant approximation" of the Oseen-Zocher-Frank free energy ([21]). 
The variational derivative 4^- is then 

81 

— = *d*da = Aa, (4.118) 
da 

i.e., simply the Laplacian of a. 
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The Continuous Kelvin-Noether Theorem. Let C be a direct sum consisting of the space of loops 
in M together with the space of scalar functions on M. Let G = Diff vo i(M) (S) !F{M, S 1 ) act on C via the 
pullback on each factor: (7, (j>) ■ (<p, 0) — (93*7, <p*(f)) for any loop 7 in M, any scalar function (f> £ F{M), any 
diffcomorphism tp £ Diff vo i(M), and any 9 £ T(M, S 1 ). Let K, : C x T(M, S 1 ) — > g** be the equivariant map 
defined by 

(/C(( 7 ,0),a),(wU)> = + <</>,<£), (4.119) 

where w is a one-form and ip £ F(M). The Kelvin-Noether theorem then gives the statement 

vV>+ I = f :-(t)Mt) + (¥-(tU(t)) (4-120) 




7 (t) J J 1 (t) \& a 

where v(f) is the fluid velocity field, uj(t) is the local angular velocity field, 7(f) is a loop advected with the 
flow, and (j>{t) = ip(t)*(f)(0) is a scalar function advected with the flow. 

To obtain a more specific conservation law, employ the free energy function (4.117) so that = Act 
and impose the boundary condition da\oM = 0. Let 0(0) = 1 and take the limit as the loop 7 contracts to 
a point. Using Stoke's theorem, we obtain 

d 
dt 



u(t)dx = / Aa{t)dx = 0. (4.121) 

M JM 



In other words, the total angular momentum due to micromotions is preserved along solutions to the gov- 
erning equations for nematic liquid crystal flow. (This is easier to see directly from (4.113), although the 
Kelvin-Noether approach illuminates the connection between this conservation law and the symmetry of 
the nematic liquid crystal Lagrangian under the action of Diff vo i(Af) on the second factor F{M, S 1 ) of 
G = DiS vol (M)®J r (M, S 1 ).) 



Spatial Discretization. A spatial discretization of nematic liquid crystal flow is obtained by replacing 
the continuous configuration space with 

G = 2?(M)©f2°(M, S 1 ), (4.122) 

where f2^(M, S 1 ) denotes the space of S' 1 -valued discrete zero- forms on the mesh M. Elements of f2^(M, S 1 ) 
are merely vectors in (5' 1 ) Ar , where N is the number of cells in the mesh. In the notation of Table 2, the 
group product is given by 

(<?i,0i)(<?2,0 2 ) = (<?i Q2,q* 2 9 1 + 9 2 ) tor(q 1 ,9 1 ),(q 2 ,9 2 )£G. (4.123) 

The Lie algebra of G is 

g = 3(M)©ftg(M). (4.124) 

We identify the space dual to g with 

0* = n 1 d (M)/dn° d (M) x n° d (M) (4.125) 

through the pairing 

((C b , tt), (£, V^)) = (C\ B) + (tt, V) (4.126) 
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for C b € n 1 d (M)/dn Q d {M), BeD(M), and tt,^ e C^(M), where 

(C b ,B) = Tr(C7 tT n5) (4.127) 

and 

(TT,1p) =tt t Qi(}. (4.128) 

The bracket on g is computed to be 

ad (AlW) (fl,V0 = (-^B,-Bw + A^), (4.129) 

with dual 

ad^fC*,*-) = (iUC b +skew(w7r T ),-A7r). (4.130) 

The space of advected parameters is f2°(M, S 1 ), whose elements a are acted upon from the right by 
elements (q, 9) of G in the following manner: 

a-(q,6) = q*a-6. (4.131) 

Note once again that f29(M, S 1 ) is not a vector space, but we will regard its members as real- valued 
discrete zero-forms on M when necessary in order to circumvent the need to introduce new notation and 
terminology. 

The infinitesimal action of g on f2^(M, S 1 ) is given by 

a- (A,u) = -Aa-u for (A,ui) e 0(M)(s)fi°(M). (4.132) 

The diamond operation (3.13) is computed as follows: For any a, j3 € f2°(M, S 1 ) and (A,to) € g, 

(a-(A,w),P)=-(Aa,0)-(u,0) 
= -a T A T Vip- (p,w) 
= -Tr(A T n(3a T ) - (p,u) 
= Tr(nA(3a T )- (fi,u) 
= -(skew(f3a T ),A)-((3,u J ) 

= -((skew(/?a T ),/3),(A^)), (4.133) 

implying that 

poa = (skew(^a T ),/3). (4.134) 
The spatially discretized nematic liquid crystal Lagrangian t : [D(M) ® Q^(M)] x &%{M, S 1 ) -> K is 

*(y,w,a) = + ^<a;,a;) - F d (a), (4.135) 

where : J7°(M, 5 1 ) — >■ R is a discrete approximation to the volume-integrated free energy F. Here, 
Y G 0(M) and w € il^(M) are the discrete counterparts to the continuous velocity field v <G Xdiv(M) and 
the continuous angular velocity field oj € J-(M), respectively. 
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The discrete-space, continuous-time Euler-Poincare equations (3.11) with the advected parameter a read 

~ + £ Y Y b = skew(^a T ) (4.136) 
ot oa 

da 

— -Ya = u, (4.138) 

where the term skew(wcj T ) in (4.136) vanishes by the symmetry of luuj t . 

If (q, 6) £ G is used to denote the configuration of the fluid, then the evolution of a is given explicitly by 

a = a ■ (<?, 6)- 1 = a ■ (g _1 , -q*0) = q*(a Q + 6). (4.139) 

A discretization of the free energy function (4.117) and its variational derivative is obtained by employing 
the tools of classical discrete exterior calculus. (See [18] for details.) Define 

F d {a)= l -\\da\\l, (4.140) 

where d denotes the classical discrete exterior derivative that takes real-valued functions on fc-simplices to 
real- valued functions on (k + l)-simplices, and || • ||2 is the discrete L 2 norm of discrete forms introduced by 
Desbrun and co-authors [181. The variational derivative p- is then 

L J oa 



— = *d*da = Aa, (4.141) 
oa 



the discrete Laplacian of a [18]. 



Temporal Discretization. To perform a temporal discretization of (4.136-4.138), we require a group dif- 
ference map t : g — > G to approximate the exponential. For the semidirect product G = T>(M) (S) £^(^j 
the exponential turns out to be a nontrivial extension of the exponential on 2?(M). The following lemma, 
whose proof is presented in Appendix A. 2, gives a formula for the exponential on G. 

Lemma 4.1. The exponential map exp : g — > G for the group G = 2?(M) (S) J7°(M, S 1 ) is given by 

exp(t(A,w)) - (e tA ,A-\l -er tA )Lj), (4.142) 

where t € M, A € 5(M), uj € 57° (M), e tA is the usual matrix exponential, and j4 _1 (J — e~ tA ) is to be regarded 
as a power series 

A-\I - e- tA ) = tl - l -± + t^- - . . . (4.143) 
2 

(so that it is defined even for A not invertible) . 

In the following three lemmas, we derive a structure-preserving approximant r to the exponential on G 
and compute its inverse right-trivialized tangent dr _1 , as well as the dual of dr" 1 . Proofs of these lemmas 
are presented in Appendix A. 2. 
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Lemma 4.2. The map t : g — > G given by 

t(A, u) = (cay(A), csy{-A/2)u) (4.144) 

is a group difference map for the group G = 2?(M) (S) f2"(M, S 1 ). That is, t is a local approximant to exp 
with t(0, 0) = (J, 0) ; and r satisfies 

t^lj)- 1 =t(-A,-uj) V(i,w)eg. (4.145) 

Lemma 4.3. The inverse right-trivialized tangent dr^ 1 :gxg->g of the map (4-144) * s given by 

dT (A,^)( B > ^) = [dc&y^B, cay(-A/2)V + ^cfcay A/2 (dcay A 1 B)w^ . (4.146) 

Lemma 4.4. For the map r given by (4-144): the dual of the operator dr~ l is given by 

>*) = (Wy^)*^ + i(dcay A 1 )*(dcay A/2 )*skew(7r W T ), cay(^/2)7r) . (4.147) 

Note that if -k is a scalar multiple of lo, then this reduces to 

(dr^TiC^n) = ((dcay^ 1 )*^, cay(^/2)7r) . (4.148) 

We are now equipped to write down the discrete-space, discrete-time Eulcr-Poincare equations (3.35-3.36) 
for nematic liquid crystal flow: Using the fact that £ = w, we use (4.148) obtain 

(dcay:jt n )*n b = (d^YhY k _J* Y k-i + hskew (jc7 k a *) (4 ' M9) 

S£ 

c&y(-hY k /2)uj k = cay(Wr fc _i/2)w fc _i + h— (4.150) 

da k 

a k = cay(/iYfc_i)(a fc _i + hcay(-hY k - 1 /2)uj k - 1 ). (4.151) 

Upon rearrangement, these become 

n - YLi , ^n-Ai + £y,n , h (Yk _ lY >_ inYk _ in -i _ r^nnn- 1 ) = skew ( (4.152) 

(4.153) 



w fc = cay (/life /2) 



cay(ftYfe_i/2)a; fe _i + /i 



a fe = cay(/iYfc_i)a fe _i + hcay(hY k _ 1 /2)uj k - 1 . (4.154) 

The Discrete Kelvin-Noether Theorem. Take C to be the Lie algebra of G = D(M) © ^"(M, S* 1 ) and 
let G act on C = &(M) © f2°(M) via the adjoint representation. Let JC : Cx Cl%M, S 1 ) -> Q be the equivariant 
map ((r,^),cn) i — ^ (T,ip). The discrete Kelvin-Noether theorem then says that the quantity 



M, 



* = ((^-in,^))*^'^)' ( r ^)) (4-155) 



4.2 Discretization of Continuum Theories 



40 



satisfies 



ith 



(r fc ,^fe) = Ad(j yfc:efc )(r ,V'o) 

= ((y k )*r ,y k (-r Q d k +-0 O )) . (4.157) 

To obtain a more specific conservation law, consider the group difference map (4.144) together with the 
free energy function (4.140). Choose To = and ipo = 1 := (1, 1, • • • , 1) T so that = and ip k = j/fel = 1 
by the stochasticity of y k ■ It then follows that 

(cay(-hY k /2)u k , 1) - (cay(-/iy fe _ 1 /2)w fe _ 1 , 1) = h(Aa k , 1) 

= 0, (4.158) 

where the last equality follows from an application of discrete Stoke's theorem to the discrete zero-form 
Aafc. Indeed, pairing such a form with the quantity 1 returns, in the framework of classical discrete exterior 
calculus, the integral of d*da k over the domain. Recasting this area integral as a line integral of *da k along 
the boundary of the domain, this quantity vanishes provided the appropriate boundary condition da k = 
is enforced along the domain's boundary. 

Using the stochasticity of cay(— hY/2), (4.158) can be rewritten as 

(w fc) l) = <w fc _ 1 ,l> (4.159) 

In other words, the total angular momentum due to micromotions is preserved exactly along solutions to the 
discrete-space, discrete-time nematic liquid crystal equations (4.152-4.154). 

Remark. Alternatively, just as we observed in the continuous case, the angular momentum conservation 
law (4.158) can be proven directly by pairing the left- and right-hand sides of (4.150) with the vector 1. 
The Kelvin-Noether approach connects this conservation law with the symmetry of the discrete Lagrangian 
under the action of G on the second factor fl d (M, S 1 ) of G = D(M) ® fi3(M, S 1 ). 



4.2.4 Microstretch Fluids 

A microstretch fluid is a complex fluid whose constituent particles undergo three types of motion: translation, 
rotation, and stretch. The variables of interest include the fluid velocity v 6 X d i v (M'), the local angular 
velocity field v € J"(M', M 3 ), the director n e J"(M',R 3 ), and the local stretch rate v € F(M',M). The 
first three of these variables have the same physical interpretation as with nematic liquid crystals. The 
new variable v§ measures the local rate of expansion of fluid particles, taking positive values at locations of 
expansion and negative values at locations of contraction. 

A final variable required for a complete description of microstretch fluid flow is the microinertia tensor 
i G J-(M' , Sym(3) + ), which assigns a symmetric positive-definite inertia tensor to each point in the domain. 
(The classic director approach to nematic liquid crystal flow presented in the previous section implicitly 
assumes i to be uniformly the identity matrix, which fortunately holds for all time if it holds at t = in the 
absence of stretching.) 
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In practice, it is common to work with the modified microinertia tensor j and its trace jo, defined as 
follows: 



3 
jo 



Tr(i)7 — i 
Tr(i) 



(4.160) 
(4.161) 



The Continuous Equations of Motion. In terms of the variables defined above, the equations of motion 
for three-dimensional, incompressible, homogeneous microstretch fluid flow are given by 



dv 



+ (v • V)v 



dF 
dn i 



Vn 



jo f dvy 

2 \ dt 
dv 
~dt 



+ v • Vt'o - Vq j + (jv) ■ v = -h • n 

(v • V)i^ — 2v jv — (jv) x v = h x n 

- (Vn)v + n x v — Vqii = 
(v • V)j + 2v j + [j, v} = 



dn 

at 



V-v = 0, 



where p is a scalar pressure field, 



dF 
dn 



di 



dF 
dn A 



(4.162) 
(4.163) 
(4.164) 

(4.165) 

(4.166) 
(4.167) 

(4.168) 



F(n, Vn) is the free energy, and v 6 J r (Af',so(3)) is the skew-symmetric matrix- valued map related to 
v € J"(M',]R 3 ) via the hat isomorphism (4.2). 

The equations of motion for microstretch fluid flow have an interpretation similar that of the governing 
equations for nematic liquid crystal flow: Equation (4.162) is Euler's fluid equation with additional terms 
arising from micromotions; equations (4.163-4.166) describe the advection of the local stretch rate Uq, the 
local angular velocity field v, the director field n, and the modified microinertia tensor j, respectively; and 
equation (4.167) expresses the incompressibility of the flow. 

Once again, we will specialize to the case of two-dimensional microstretch fluid flow for simplicity. Assume 
that the motion is restricted to the x-y plane so that the director n is everywhere orthogonal to the z-axis 
and the modified microinertia tensor j has the form 



(4.169) 





(in 


J12 


o\ 


3 = 


hi 


J22 







V° 





333/ 
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Define the variables (ui,R) £ J"(M,M © K), (a, A) £ J"(M, S 1 © R), (j r> j s ) £ JF(M, R®R) via 

i/ = (0,0,w) T (4.170) 

= i? (4.171) 

n = e A (cosa,sina,0) T (4.172) 

^io = (4.173) 

J 33 = e^. (4.174) 

The motivation behind this choice of variables will be made clear in the next section. 

The governing equations for two-dimensional, incompressible, homogeneous microstretch fluid flow are 
then given by 

% + ■ v >* = -* O v ° - a < (19 va - v " <« 75 > 

W +V ' V ^^fcj-^ (4 - 176) 

<7Q + v • VQ = -2z r - 2* s + 9, - ^ (4.177) 



v 



dt ^ ' 51 \d\,iJ d\ 

da 
~dt 
dX 
dt 

djr 

dt ' 
dj s 
dt ' 



vVa = w (4.178) 
• VA = R (4.179) 



v • Vj r = -2R (4.180) 
v • Vj, = -2R, (4.181) 



where 



7T = e^uj (4.182) 
Q = e j °R (4.183) 



= Vw 2 



(4.184) 



= ^e Js i? 2 . (4.185) 



Geometric Formulation. In yet another remarkable fashion, the governing equations for microstretch 
fluid flow can be cast as Euler-Poincare equations on a certain Lie group with respect to a natural Lagrangian. 
In this case, as Gay-Bahnaz & Ratiu [21] show, the appropriate configuration group (on a three dimensional 
domain M') is the semidirect product Diff vo i(M') (S) J-(M', CSO(3)), where CSO(3) is the conformal special 
orthogonal group: the group of (positive) scalar multiples of rotation matrices. 

In two dimensions, the configuration space for microstretch fluid flow is the semidirect product 

G = Diff vol (M) ® F(M, S 1 © R), (4.186) 
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where J"(M, S 1 © R) is the space of maps from the domain to the direct sum S 1 ® R. The group product is 
given by 

{<pi,{9 1 ,r 1 ))(ip 2 ,{02,r 2 )) = (yi 0(^2,(^1 +6 , 2,<P^i + r 2 )) for (6>i, n)), (y> 2 , (6> 2 , r 2 )) G G. (4.187) 
We regard G as the two-dimensional realization of the group 

G' = Diff vo i(M') ®T{M', 050(3)) (4.188) 

with group product given by 

(Vi> Xi)(<P2, X2) = (<Pi o ¥> 2 , (V2XOX2), (4.189) 

where M' is an open subset of R 3 containing an embedding of M in the x-y plane and any \ G T[M' ', CSO(S)) 
is related to a quantity (0, r) G J"(M, S" 1 ® R) via 



(4.190) 





^ cosO 


— sin 
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V 



The Lie algebra of G is 

= Xdiv(M) © J"(M, 101). 

We identify the space dual to Q with 

q* = fl 1 (M)/dfl (M) x J"(M,R®R) 

through the pairing 

((w b , (tt, 5)), (v, R))) = (w b , v) + (tt, V) + (5, i?> 
for w b e n 1 (M)/dn°(M), v G Xdiv(M), and (tt,S), (tp,R) G J"(M,R® R), where 

(w b ,v) = / w b (v)dx 
(ir,ip) = / 7r(x)^>(x) dx 
(5, i2) = f S(x)R(x) dx. 

The bracket on g is computed to be 

ad( u ,(«,o))(v, (V>, #)) = (- AiV, (£ v u - £ u 4>, £ V Q - £uR)), 

with dual 

ad (u,(^,Q))( wb > s )) = ( £ uW k + nduj + SdQ, (£ u ir, £ U S)). 

The advected parameters are (a, A) G J*(M, S 1 ® R) and (>, j s ) G J*(M,R® R) 
from the right by elements (ip, (9, r)) of G via 

(a,A)-(<^ (6,r)) = (<p*a-6,v*\-r) 
(jr,j s ) ■ {ip, (6, r)) = (<£*> + 2r, ^*j s + 2r). 



(4.191) 
(4.192) 
(4.193) 

(4.194) 
(4.195) 
(4.196) 

(4.197) 
(4.198) 

which are acted upon 

(4.199) 
(4.200) 
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The choice of such a set of advected parameters and the corresponding group actions (4.199-4.200) is mo- 
tivated by the ambient geometric structure of three-dimensional microstrctch fluid flow. In three dimensions, 
the advected parameters are the director n € T(M' , R 3 ) and the microinertia tensor i € T[M' , Sym(3) + ), 
related to (a, A) and (j r ,js) via (4.172-4.174). These parameters are acted upon from the right by elements 
(<p,X) of BiS vo i(M')®T(M',CSO{3)) via 

n-(<P,x)=X~ 1 {<P**>) (4.201) 
i-(¥>.x)=X 3 Vi)x- (4-202) 

One checks that for x of the form (4.190), the three-dimensional actions (4.201-4.202) on n and i induce 
the two-dimensional actions (4.199-4.200) on (9,r) and (j r ,js)- One furthermore checks that, in a rather 
convenient manner, the quantities j r and j s are sufficient to reconstruct the products ju and jo^o appearing 
in (4.162-4.167), provided the modified microinertia tensor j = Tr(z)J — i has the form (4.169). 

Note that the exponentials appearing in relations (4.172-4.174) play an essential role in ensuring the 
feasibility of an eventual discretization of the configuration space and its action on the advected variables. 
They endow the space of advected parameters with an additive, rather than multiplicative, structure, thereby 
guaranteeing the sensibleness of replacing continuous pullbacks of scalar fields by matrix- vector products for 
the purposes of discretization. 

Having made note of these observations, let us proceed with a derivation of the Euler-Poincare equations 
associated with the group G given by (4.186). The infinitesimal actions of q = Xdi v (M)(§)J r (M, K © M) on 
the advected variables (a, A) € J"(M, S 1 © R) and (j r ,j s ) G J"(M, R © R) are given by 

(a, A) • (u, (cj, R)) = (£ u a - u, £ U X - R) (4.203) 
(j r ,j s ) ■ (u, (w, R)) - (£ u j r + 2R, £ u j s + 2R). (4.204) 

These give rise to the following diamond operations, which are easily derived through analogy with (4.109): 

(/3, a) o (a, A) = (-fida - od\, (/3, a)) (4.205) 
(i r ,i s ) o [jr,j s ) = {~i r dj r - i s dj s , (0,-2i r - 2i s )). (4.206) 

The Lagrangian £ : [X div (M)®J r (M, R © R)] x {F(M, S 1 © R) © J"(M, R © R)] -)• R for two-dimensional 
microstretch fluid flow is given by 

£(v 1 u Jl R,a,\ 1 j r J s ) = l(v l ',v) + l(e^L},u J ) + l(eJ°R,R)- [ F(a, da, A, d\)dx, (4.207) 
* * * Jm 

where F is the microstretch free energy function. The first three terms in this Lagrangian correspond to, 
respectively, the kinetic energy due to translational, rotational, and stretching motion. 
Make the following definitions for the variational derivatives of £: 

- §L 

Q = - 

W SR 

51 





(4.208) 


e 3 'R 


(4.209) 


1 , o 
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The Euler-Poincare equations (3.11) associated with the Lagrangian (4.207) then read 

^ + £ v v b =-~da-~d\-dp (4.212) 
ot oa oa 

% + £ v n = £ (4.213) 
Ot da 

d § r + £ v Q = -2i r - 2i s + ~ (4.214) 

Ot OA 

da 

— + i? v a = w (4.215) 
ot 

d\ 

— + £ V \ = R (4.216) 
or 

% + £vir = -2i? (4.217) 
ot 

% + £ v i s = -2R, (4.218) 

where four of the terms in (4.212) have conspired to produce a full differential, which we have absorbed into 
the pressure differential dp. 

A convenient choice for the free energy F is 

F(a 7 da,X,d\) = ^\\da\\ 2 + ^A 2 , (4.219) 

which corresponds to a modification of (4.117), obtained through the addition of a potential energy due 
to internal stretching in accordance with a leading order approximation to Hooke's Law. (Recall that the 
length of the director n is related to A via ||n|| = e A .) The variational derivatives and M- are then 

Si 

— =*d*da = Aa (4.220) 
oa 

S£ 

- = -A. (4.221) 

The Continuous Kelvin-Noether Theorem. In analogy with nematic liquid crystal flow, the Kelvin- 
Noether theorem applied to microstretch fluid flow with the free energy (4.219) and the boundary condition 
da\dM = gives, as one corollary, the conservation law 



d 
dt 



/ n(t)dx = 0. (4.222) 

J M 



That is, the total angular momentum due to micromotions is preserved along solutions to the governing 
equations for microstretch fluid flow. 

Spatial Discretization. A spatial discretization of two-dimensional microstretch fluid flow is obtained by 
replacing the continuous configuration space with 

G = 2?(M)(s)^(M,S' 1 ®M), (4.223) 
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where 05(M, S 1 © M.) denotes the space of (S 1 © R)-valued discrete zero-forms on the mesh M. For a mesh 
with N cells, elements of S1°(M, 5 1 © R) are merely pairs (6,r) of real-valued N X 1 vectors and r, with 
the entries of 9 taken modulo 2ir. The group product is given by 

(q h (fli,ri))(ga, (0 2 , r 2 )) = (g iga , % + 9 2 , q^n + r 2 )) for (0i,n)), (©, (0 2 ,r 2 )) € G. (4.224) 

The Lie algebra of G is 

fl = 3(M)(S)f^( M : K © K )- (4.225) 

We identify the space dual to g with 

g* = n 1 d (M)/dn° d (M) x fi°(M,R©R) (4.226) 

through the pairing 

((C b , (tt, S)), (B, (V, fl))) = (C b ,B) + (vr, + (5, i£) (4.227) 

for C b G nJ(M)/dng(M), B G D(M), and 7r,ip,S,R G ^( M ), where 

(C b ,B) = Tr(C bT flB) (4.228) 

=n T flxjj (4.229) 

(5, fl) = S* T fti?. (4.230) 

The bracket on g is computed to be 

ad (X ,( u ,Q))(B,(^,iJ)) - (-.^(-Bw + A^-BQ + AR)), (4.231) 

with dual 

ad( A , (w> Q))(C b , fa S)) = (£ A C b + skew(w^ T ) + skew(QS T ), {-Air, -AS)). (4.232) 

The advected parameters are (a, A) G Cl d (M, S 1 ©R) and (>, J s ) G fi°(M, R © R), which are acted upon 
from the right by elements (q, (9, r)) of G via 

(a, A) ■ (q, (9, r)) = (cT 1 ** - 0, ^T 1 A - r) (4.233) 
(j r , j s ) ■ (g, (0,r)) = (q-'jr + 2r,q- 1 j s + 2r). (4.234) 

The induced infinitesimal actions of g = D(M) ® fi°(M, R © R) on the advected variables (a, A) G 
fi£(M, 5,1 © K ) and (>,j s ) G ft°(M,R © R) are given by 

(a, A) ■ (A, (w, i?)) = (-Aa - w, -AA - J?) (4.235) 
C? r> j s ) ■ (A, {cj, R)) = (-Aj P + 2iZ, -Aj s + 2R). (4.236) 

The astute reader will have recognized the pattern relating continuous and discrete diamond operations 
on scalar fields discussed in Section 4.2.3, thereby permitting an easy calculation of the discrete diamond 
operations via inspection of their continuous counterparts (4.205-4.206): 

09, a) o (a, A) = (skew(/3a T ) + skew(crA T ), (J3, a)) (4.237) 
(i r ,i s ) o tir,j s ) = (skew(i r .^) + skew(i s jj), (0, -2i r - 2i s )). (4.238) 
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The spatially discretized Lagrangian I : [o(M) ® fl° d (M, R © R)] x [Q d (M, S 1 © R) © ft°(M, R © R)] -)• R 
for two-dimensional microstrctch fluid flow is given by 

£(Y,u,R,a, \,j r ,js) = l -(Y\Y) + ^Qu,u) + l -{e^QR,R) - F d (a,X), (4.239) 

where Fd '■ ^l d (M, S 1 ©R) — > R is a discrete approximation to the volume-integrated microstretch free energy. 
The notation here refers to the Hadamard (entry-wise) product of vectors. Likewise, the exponentials 
appearing here are entry-wise exponentials. 

In analogy with the continuous definitions from before, make the following definitions for the variational 
derivatives of £: 



51 
51 



e jr ui (4.240) 



Q= — =e^QR (4.241) 
OR 

se i ■ 

i r = — = ~e 3r 0w0 w (4.242) 

Ojr 2 

51 I ■ 

i s = — = -e J ° &R&R. (4.243) 
ojs 2 

The discrete-space, continuous-time Euler-Poincare equations (3.11-3.12) associated with the Lagrangian 
defined in equation (4.239) then read 

+ £yY b = skew(7ra; T ) + skew(Qi? T ) + skew(vj'J) 



at 



+ skew(i s jj) -I- skew (J^^j + skew (jf^^^) (4.244) 

%-Y,= f- (4.245) 
ot da 

^-YQ = -2%r - 2i s + ^ (4.246) 



da 
~dt 
dX 
dt 

djr 

at ' 
dj s 

dt ' 



- Ya = cu (4.247) 

-Y\ = R (4.248) 

Yj r = -IR (4.249) 

Yj s = -2R. (4.250) 



Note the presence of the four extra terms in (4.244) relative to (4.212); these terms only conspire to produce 
a full discrete differential in the continuous limit. 

A discretization of the free energy (4.219) is given by 

F d (a,A) = i||da||l + i{A J A). (4.251) 
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The variational derivatives ¥- and |i are then 

da d\ 

Si 

— =*d*da = Aa (4.252) 
oa 

fx = x ^ 

Temporal Discretization. A temporal discretization of (4.244-4.250) is obtained most easily by choosing 

t(A,u) = {csy{A),u), (4.254) 

which is not a group difference map in the sense of Definition (3.4) (since t(A,lu)^ 1 ^ t{—A, -co)), but it 
nonetheless maps q into G and locally approximates the exponential. While this choice provides a lower 
order approximation to exp than the genuine group difference map (4.144), it has the appealing feature that 
its inverse right-trivialized tangent has a simpler form than that shown in (4.146). 

Note that in general, the use of a map r with t(£) -1 ^ r(— £) for £ 6 g leads to the following generalization 
of the (right-right) discrete Euler-Poincare equations with an advected parameter (3.34-3.36): 

g k+ i = r(h^ k )g k (4.255) 

a k+1 = a fe T(/i£ fc )- 1 (4.256) 

(Dr-\r(h^))-TL T{hik) y— = (Dt^t^-i)) • TR^^)* — + h— o a k , (4.257) 

ot, k oa k 

which are easily derived through an argument that parallels the proof of Theorem 3.7. 
For the map (4.254), one computes that for (A,uj) <E Q, (C ,tt) € Q* , 

(Dt-\t(A,lo)))-TL t(a ^ )) )*{C\k)= ((dcay:\)*C b -skew(7r W T ),^) (4.258) 

(D T - 1 (T(A,uj)))-TR T(AtUj)) )*(C l ',ir) = ((dcay^ 1 )*^, cay(A)7r) (4.259) 

Substituting into (4.256-4.257) gives the discrete-space, discrete-time Euler-Poincare equations for mi- 
crostretch fluid flow: 

= skew(7r fc Wfc) + skew(Q k R{) + skew(i rifc j^ fc ) + skew(i S!fe j^ fc ) 

+ skew (^«fc J + skew (J^^j ( 4 - 260 ) 
81 

n k = cay(/iF fe _i W_i + h— (4.261) 

da k 

Si 

Q k = cay(ft,Yfe_i)Q fc _i - 2hi r . k - 2hi s t + /i^r- (4.262) 

a fe = cay(ft,r fe _i)(o! fe _i + hu k -i)- (4.263) 

A fc = cay(/ l F fe _ 1 )(A fc _ 1 + hR k _ x ). (4.264) 

ir.fc = cay(Wfc_i)(j r , fe _i - 2hR k -i). (4.265) 

j s , fc = cay(/iYfc_i)(j Sifc _i - 2hR k -i). (4.266) 



5 Update Equations on a Cartesian Grid 



49 



The Discrete Kelvin-Noether Theorem. The Kelvin-Noether theorem applied to microstretch fluid 
flow with the discrete free energy (4.251) and the boundary condition da\dM = gives, as one corollary, the 
conservation law 

(7r fc ,l) = <7r fc _i,l). (4.267) 

That is, the total angular momentum due to micromotions is preserved exactly along solutions to the discrete- 
space, discrete-time microstretch continua equations (4.260-4.266). 

5 Update Equations on a Cartesian Grid 

We now specialize to the case in which the mesh M is a two-dimensional cartesian grid and derive the 
cartesian realizations of the discrete update equations (4.51), (4.82-4.83), (4.152-4.154), and (4.260-4.266) 
for ideal fluid flow, MHD, nematic liquid crystal flow, and microstretch fluid flow, respectively. While the 
matrix forms of the update equations are readily implementable (on regular or irregular meshes), the reader 
will find that upon specializing to a cartesian grid and analyzing the matrix equations in an entry-wise 
fashion, easily implementable finite difference schemes emerge with a more familiar form. 

In fact, a comparison with the literature will reveal that our MHD scheme, in particular, employs a spatial 
discretization (although not its temporal discretization) that is identical to that of the MAC-Yee scheme on 
a two-dimensional cartesian grid [32]. As a method that integrates two popular structure-preserving finite 
difference schemes the MAC scheme from fluid dynamics and the staggered mesh-based Yee scheme from 
electrodynamics - the MAC-Yee scheme debuted in 2001 as a promising integrator for ideal incompressible 
magnetohydrodynamics that maintains the V • B = constraint exactly and respects energy and cross- 
hclicity conservation in the following sense: the spatially discretized, continuous-time equations presented 
by Liu and Wang [32] conserve a discrete energy and a discrete cross-helicity. 

Our MHD integrator surpasses these milestones in several key respects. Firstly, while our spatial dis- 
cretization coincides with that of the MAC-Yee scheme, our temporal discretization differs markedly. As a 
consequence of our careful choice of temporal discretization: 

1. Our integrator is variational, suggesting that our integrator should exhibit good energy behavior even 
after temporal discretization, not merely in the context of the spatially discretized, continuous-time 
equations. 

2. Through the discrete Kelvin-Noether theorem, our integrator preserves a discrete cross-helicity exactly 
after temporal discretization. 

Finally, our derivation of an MHD integrator was geometric, providing physical insight into the success of 
staggered grid schemes and generalizing the MAC-Yee scheme to unstructured meshes. 

The Discrete DifFeomorphism Group on Cartesian Grids. In the case of a two-dimensional cartesian 
grid M with uniform spacing e, the diagonal matrix of cell volumes Q coincides with a multiple of the identity 
matrix: 

n = e 2 I. (5.1) 
Consequently, the group P(M) reduces to the group of orthogonal, signed stochastic matrices: 



X>(M) = {qe GL(N)+ | Qij = 1 Vt, q T q = I}, 

3 



(5.2) 



5 Update Equations on a Cartesian Grid 



50 



and its Lie algebra Z)(M) reduces to the space of antisymmetric, row-null matrices: 



0(M) S 1(N) I A H = Vi,A T + A = 0}. 



(5.3) 



3 



Matrices A in the nonholonomically constrained space S C 0(M) that approximate vector fields v = 
(u, v) £ Xdi v (M), McR 2 , have nonzero entries given by 



where D mn denotes the edge shared by cells C rn and C n and x m „ £ M is the position of its midpoint. 

The constraints (5.3) appearing in the definition of 0(M) have the following implication for matrices 
A G Z)(M) approximating continuous vector fields: the fluxes of the velocity field across the edges bounding 
any given cell C n must sum to zero. This constitutes the discrete analogue of the constraint V • v on the 
velocity field in incompressible fluid flow, and we will record this constraint in every finite-difference scheme 
presented below for completeness. 

We should emphasize, on the other hand, that the MHD cartesian update equations (5.41-5.42) (to be 
derived in this section) automatically preserve the divergence-freeness of the magnetic field at the discrete 
level since the ambient matrix update equation (4.81) is an equation on the Lie algebra of antisymmetric, 
row-null matrices, a space that is closed under addition and commutation. This is consistent with the fact 
that in the continuous world, V • B holds automatically if it holds at t = 0, whereas V ■ v = can be viewed 
as a constraint that uniquely determines the pressure field p. 

Nonholonomic Constraints and Cubic Terms. As discussed in Section 2.3, we have imposed nonholo- 
nomic constraints throughout Section 4.2 to guarantee the desired sparsity of the matrices used to represent 
vector fields and one-forms. When dealing with the resulting weak equalities that appear in (4.51), (4.82- 
4.83), (4.152), and (4.260), the following two procedures will be employed, in accordance with the remarks 
in Section 2.3: 

1. When working with Lie derivatives £aC^ of discrete one-forms C b € fl d (M)/dfl d (M.), we concern 
ourselves only with those components (£AC^) mn of £aC^ corresponding to adjacent cells C m and C n . 

2. When working with discrete vector fields in D(M), we apply the sparsity operator (2.36) to any com- 
mutators that arise. 

Lastly, for simplicity, we shall drop any cubic terms (matrix products involving three elements of 0(M) U 
n 1 d (M)/dn° d {M)) appearing in the update equations (4.51), (4.82-4.83), (4.152-4.154), and (4.260-4.266). 
This practice is advocated in the context of rigid body simulations by Kobilarov and co-authors [29], who 
point out that the dropping of such terms does not reduce the second-order of accuracy of the discrete Euler- 
Poincare equations arising from the Cayley group difference map. Note also that the discrete Kelvin-Noether 
theorem still holds exactly if dcay^ 1 (Z) and its dual (dcay^ 1 )* X b are replaced by 




(5.4) 




(5.5) 



5.1 Frequently Encountered Quantities and their Cartesian Realizations 



51 



and 

1 



I--£ Y )X\ (5.6) 
respectively - which is equivalent to dropping the cubic terms. (Compare with (4.48) and (4.50).) 

Notation. In what follows, we parameterize the vertices of our coordinate grid with integer pairs G 
Z x Z. Spatial locations of quantities will be encoded with superscript notation. Thus, vertex-based quantities 
are denoted , cell-centered quantities are denoted gr i + 1 / 2 ^*+ 1 / 2 J and edge-based quantities are denoted 
qi+i/2J anc [ q,jj+i/2_ a s shorthand for averages between pairs of adjacent cell-based quantities, we shall use 
the following intuitive notation: 

q i+l/2,j ._ 1 ^1+1/2^-1/2 + q i+l/2,j + l/2^ £g _^ 

and similarly for averages between pairs of adjacent edge-based quantities and pairs of adjacent vertex-based 
quantities. In some instances, products of compatible quantities will be combined if there is no danger of 
ambiguity, e.g. 

(grji+i/aj := g ;+i/2,v+i/2j (5 8) 

In a slight abuse of notation, subscripts will be used for both matrix/vector indexing and temporal indexing. 
While context should alleviate most ambiguities, we shall nonetheless reserve the subscript k for temporal 
indexing; any other subscripted letter shall indicate componentwise indexing of matrices and vectors. 

Finally, the discrete curls of vector fields u of the form u = (u, v) shall play an important role; these 
vertex-based quantities are defined as follows: 

u iJ-l/2 + v i+l/2,j _ u i,j+l/2 _ v i-l/2,j 

uj hJ (u,v) = . (5.9) 

Fig. 5.1 provides a useful reference for the subsequent discussion. 

5.1 Frequently Encountered Quantities and their Cartesian Realizations 

Below we summarize the cartesian realizations of some frequently encountered quantities arising in discretiza- 
tions of flows on the diffeomorphism group. Only results are stated here; derivations of these results are 
presented in Appendix A. 3. 

5.1.1 The Lie Derivative of a Discrete One-Form 

We will begin by considering the cartesian realization of quantities of the form £a& with A, C € S C 0(M). 
This quantity has, not surprisingly, appeared in every continuum theory discretization presented in this 
report. It is the basic object governing flows on the discrete diffeomorphism group, serving as the discrete 
counterpart to continuous quantities of the form £ u w with u, w G Xdi v (-^)- 

Assume that the discrete vector fields 4,C £ 5 approximate continuous vector fields u = (u, v),w = 
(r, s) € %di v (M), respectively, so that the entries A mn and C mn for a pair of horizontally adjacent cells C m 
and C n centered at (i — 1/2, j + 1/2) and (i + 1/2, j + 1/2) are, in accordance with (2.25), given by 

A mn = -^ +1 ' 2 (5.10) 

C mn = -j/" J+1/ \ (5.11) 
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'Jyi'-1/2J+J 


u lJ+m 




u'- j+m 







i,j+l/2 



i+lj+l/2 



(a) 



(b) 



Figure 5.1: (a) Stencil for computations involving horizontally adjacent cells, (b) Stencil for computations 
involving vertically adjacent cells. 



Likewise, the entries A mn and C mn for a pair of vertically adjacent cells C m and C n centered at (i+l/2, j — 1/2) 
and (i + l/2,j + 1/2) are given by 



(5.12) 
(5.13) 



Fix a pair of vertically adjacent cells C m and C n . Then it can be shown (see Appendix A. 3) that, modulo 
a discrete differential, 



(£ A C%„ = -~ (w«(r, S ) f — +o,^(r,.) — 



(5.14) 



Likewise, for a pair of horizontally adjacent cells C m and C n centered at (i — 1/2, j + 1/2) and (i + 1/2, j - 
1/2), we have 



(£ A C") mn = - w w (r,s +^< J+1 (r, s ) 



(5.15) 

Notice that we have recovered, at the discrete level, the two-dimensional manifestation of the vector 
calculus identity 

£ u w b = ((V x w) x u) b , (5.16) 
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which holds for vector fields u,w defined on M 3 , modulo a full differential. 



5.1.2 The Lie Derivative of a Discrete Vector Field 

In this section, we consider the cartesian realization of quantities of the form £aB with A, B G S C 0(M). 
Quantities of this form arise in the discrete Euler-Poincare equations (4.77-4.78) for magnetohydrodynamics; 
they serve as the discrete counterparts to continuous quantities of the form <£ u v with u, v € X^wiM). 

Assume that the discrete vector fields A, B E S approximate continuous vector fields u = (u,v),v = 
(p, q) € Xdi V (M), respectively, so that the entries A mn and B mn for a pair of horizontally adjacent cells C m 
and C n centered at (i — 1/2, j + 1/2) and (i + 1/2, j + 1/2) are given by 

A mn = ~u i ^ 2 (5.17) 

B mn = -^+ 1 ' 2 , (5.18) 

Likewise, the entries A mn and B mn for a pair of vertically adjacent cells C m and C n centered at (i+1/2, j — 1/2) 
and (i + 1/2, j + 1/2) are given by 

A mn = ~v i+1 ' 2 * (5.19) 

B mn = -~q i+1 / 2 ' j . (5.20) 

Fix a pair of vertically adjacent cells C m and C n . Then, in the abbreviated notation discussed in this 
section's introduction, it can be shown (see Appendix A. 3) that 

. 1 /yi+ljqi+iJ - u ijqij //•' .•',■" ' ' />••■'<•' '\ 

(£AB) mn — —— I - - I . (5-21) 

Similarly, for a pair of horizontally adjacent cells C m and C n centered at (i— 1/2, j + 1/2) and (i + 1/2, j + 
1/2), we have 

Notice that here we have recovered, at the discrete level, the two-dimensional manifestation of the vector 
calculus identity 

£ uV = V x (v x u), (5.23) 
which holds for divergence- free vector fields u, v defined on R 3 . 

5.1.3 The Antisymmetrization of an Outer Product 

Below we consider the cartesian realization of quantities of the form skew(F_E T ) £ f7^(M)/d57°(M), where 
E and F are discrete zero-forms. As the reader of Sections 4.2.3-4.2.4 is well aware, such quantities arise 
profusely in discretizations of scalar field advection by way of the diamond operation. See for instance, 
the Euler-Poincare equations (4.136-4.138) for two-dimensional nematic liquid crystal flow and the Eulcr- 
Poincare equations (4.244-4.250) for two-dimensional microstretch fluid flow. 
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Assume that the discrete zero- forms E,F E f^(M) approximate continuous scalar fields a, P € T(M), 
respectively, so that the entries E n and F n for a cell C n centered at (i + 1/2, j + 1/2) are given by 

E n = a *+l/2 J + l/2 (5 24) 

F n = ^+V2,i+i/2, (5.25) 
Recall the correspondence between skew(_Fi5 T ) and its continuous counterpart: 

Continuous Discrete 

[-fida] [skew(F£ T )] (5 . 26) 

m m 

n l (M)/dn°(M) n 1 d (M)/dn° d (M) 

Here, we have explicitly included the brackets to emphasize the distinction between cosets and their repre- 
sentives. 

Computing the (m, n) entry of skew(FE T ) for vertically adjacent cells C m and C n centered at {i + 1/2, j — 
1/2) and (i + 1/2, j + 1/2), respectively, one finds (see Appendix A. 3) that 

e/ / fl*+l/2,J+l/2 _ oi+l/2j-l/2\ / ^+1/2^ + 1/2 _ ^+1/2,^-1/2 



2 V V 6 

(5.27) 

Similarly, for a pair of horizontally adjacent cells C m and C n centered at (i— 1/2, j + 1/2) and (i + 1/2, j + 
1/2), one obtains 

e/ , / fl*+l/2,J+l/2 _ m-l/2j + l/2\ /n,<+l/2,i+l/2 _ ^-1/2,^+1/2 

skew(Fi? T ) m „ =-£( 0^+^ £ £ _ ^+1/2 



2 v V e 

(5.28) 

Evidently, the matrix skew(FE T ) encodes the continuous one-form ^(ad/3—f3da). This agrees with (5.26), 
up to the choice of representative for the coset [— (3da\. 

5.2 Summary of the Cartesian Update Equations 

Having identified the cartesian realizations of the three most frequently encountered quantities appearing in 
the Section 4's discretizations, we now record the cartesian update equations for ideal fluid flow, MHD, liquid 
crystal flow, and microstrctch fluid flow in their most compact form, suitable for numerical implementation. 

In anticipation of the ubiquity of terms of the form (5.14-5.15) with (u, v) — (r,s), define the following 
two operators, which take a cartesian (edge-based) discrete vector field (u, v) to a cartesian (edge-based) 
discrete vector field (^ x , ^y)'- 

*-+ 1/2 M = -5 ( ( P ) +^+\u,v) ( V - P ) ) (5.29) 



2 V V 2 / V 2 

(«, «) = - (w« [u, v) 1 J + iS* 1 ** (u, v) [ - J J . (5.30) 
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Likewise, for cartesian (edge-based) discrete vector fields (u,v) and (p,q), define 

^/\u,v,p,q) = q - q + P - P (5.31) 

&+^(u,v,p, q ) = q U -^- P - (5.32) 



Note that this definition employs the abbreviated notation for averages introduced in this section's intro- 
duction by writing, for instance, u 1 ^ as shorthand for ^(u*'- 7-1 / 2 + it 1 '- 7 " 1 " 1 / 2 ). 

In addition, for cartesian (cell-based) discrete zero-forms a and /3, define the following cartesian (edge- 
based) discrete vector field (A x , A y ): 



1 / , / fli+l/2j + l/2 _ m-l/2J+l/2\ / i+l/2,j+l/2 _ _,i-l/2,j+l 
A^ + 1 /2 (/3;a) = 1 r a i^+l/2 r P _^ j _ ^ + 1/2 f « _^ 

1 / / fl*+l/2,J + l/2 fli+l/2,j-l/2\ 



(5.33) 



»+l/2,j+l/2 _ 0,1+1/2^-1/2 



(5.34) 



where once again we have employed abbreviated notation for pairwise averages between cells. 
Lastly, let A(a) denote the cartesian (cell-based) discrete Laplacian of a discrete zero-form a: 

i-l/2,j+l/2 , i+l/2,j-l/2 , i+3/2j + l/2 , i+l/2,j+3/2 _ a i+l/2,j+l/2 

A «+i/2 J+ i/2( a | = TJ± Tj" _ ™ . (5.35) 

In terms of these operators, the update equations for ideal fluid flow, MHD, liquid crystal flow, and 
microstretch fluid flow on a two-dimensional cartesian grid with spacing e and time-step h are now reviewed. 

5.2.1 Ideal, Incompressible Fluid Flow 

Variables: 



Variable 


Meaning 


Residence 


Indexing 


u 


Velocity field, horizontal component 


Vertical edges 


u iJ+l/2 


V 


Velocity field, vertical component 


Horizontal edges 




P 


Pressure 


Cell centers 


p i+l/2,j+l/2 



Update equations: 

ij+l/2 ij+1/2 lTl i,j+l/2/ \ . , T ,i,j+l/2, \ i+l/2,j+l/2 i-l/2,j+l/2 

u k ~ U k-l . *x (ttfe_i,Ufe_i) + *x' J (Ufe,Ufe) _ P k -Pk 



(5.36) 



h + 2 _ e (5 - d7) 

»+l,j+l/2 _ i,j+l/2 i+l/2,j+l _ i+l/2J 

^ ^ + ^ ^ = (5.38) 
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5.2.2 Ideal, Incompressible Magnetohydrodynamics 

Variables: 



Variable 


Meaning 


Residence 


Indexing 


u 


Velocity field, horizontal component 


Vertical edges 


u iJ+l/2 


V 


Velocity field, vertical component 


Horizontal edges 




P 


Pressure 


Cell centers 


p»+l/2,j+l/2 


B x 


Magnetic field, horizontal component 


Vertical edges 


R»,j + l/2 

J->x 


By 


Magnetic field, vertical component 


Horizontal edges 


R *+l/2,j 

Oy 



Update equations: 

# ^fc U k-1 _j_ ^x 



(itfe_i,Ufc_i) + *L J+1/2 (w fc ,w fc ) 



, T /J + l/2/ R R U* ij ' +1 / 2 m R ^ *+V2j + l/2 »-l/2,j+l/2 

x (,-Dx,fe-lj &y,k-l) + *x \-t>x,k, By.k) Pk ~ Pk 



i+l/2.j _ „i+l/2,j ^i+l/2,j 



'K-i,«fc-i) + ^ +1/2 ' J K,f fc ) 

h 2 

_ *i/ +l ^ 2 '' 7 (-8x^-15 By.k-\) + ^y^^^iBx^, By^k) P k 



»+l/2,j+l/2 _ t+l/2,j-l/2 



gi,j + l/2 _ gi,j + l/2 

B i+l/2,j _ R i+l/2,j 



Bx,k-i + B x k By k-i+By k . 
Ufc-i,Ufe_i, — : , — — I = 



i+l,j+l/2 ij+1/2 i+l/2,j+l „,»+l/2,j 



— u 



k 



+ 



k 



— V 



k 



(5.39) 

(5.40) 
(5.41) 

(5.42) 
(5.43) 



5.2.3 Nematic Liquid Crystals 

Variables: 



Variable 


Meaning 


Residence 


Indexing 


u 


Velocity field, horizontal component 


Vertical edges 


u iJ+l/2 


V 


Velocity field, vertical component 


Horizontal edges 


v i+l/2,j 


P 


Pressure 


Cell centers 


p i+l/2,j+l/2 


OJ 


Local angular velocity field 


Cell centers 


w *+l/2,j+l/2 


a 


Director orientation 


Cell centers 


a i+l/2,j+l/2 
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Update equations: 



, T ,i,j+l/2/ \ . , T ,ij + l/2, \ 



ft 2 

i+l/2J+l/2 _ i-l/2,j+l/2 

= A^ +1 / 2 (A afe , a fe ) - *± ^ (5.44) 



i+l/2J _ i+l/2,j lTr i+l/2,j 

fc fe— 1 



, Tl i+l/2,j/ \ . , T ,i+l/2,j/ \ 



i+l/2,j+l/2 i+l/2,j-l/2 

= Al+^(Aa k ,a k ) - & ^ (5.45) 



Equations (4.153 - 4.154) 1 (5.46) 



i+lj+1/2 jj+1/2 i+l/2j + l j+l/2,j 



Hi — It. V, ' " — v, , 
_fc * + _fc k = q ( 5 _ 47 ) 



5.2.4 Microstretch Fluids 

Variables: 



Variable 


Meaning 


Residence 


Indexing 


it 


Velocity field, horizontal component 


Vertical edges 


tt *j+i/a 


V 


Velocity field, vertical component 


Horizontal edges 


v i+l/2j 


P 


Pressure 


Cell centers 


p i+l/2,j + l/2 


U) 


Local angular velocity field 


Cell centers 


u i+t/2,j+l/2 


R 


Local stretch rate 


Cell centers 


R i+l/2.j + l/2 


a 


Director orientation 


Cell centers 


a i+l/2,i+l/2 


A 


Logarithm of director length 


Cell centers 


A i+l/2j + l/2 


jr 


Local rotational inertia 


Cell centers 


.i+l/2j + l/2 
Jr 


3s 


Local stretch inertia 


Cell centers 


•i+l/2,j+l/2 


IT 


*=&=e*-©w 


Cell centers 


^4+1/2,^+1/2 


Q 




Cell centers 


Qi+l/2,j+l/2 




i r = #■ = |e-?»- 0cj0w 


Cell centers 


^+1/2,^+1/2 




j s = «. = QRQR 


Cell centers 


.<+l/2,j+l/2 

OS 



x We do not recast these equations in cartesian form, as their structure does not admit such a representation. (No rearrange- 
ment of the equations can eliminate the appearance of a matrix inverse.) This is of no concern, as they are already in a form 
most suitable for efficient implementation; indeed, they are explicit updates. 
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Update equations: 



Ui 



h 2 

= A% j+1/2 (n k , u k ) + A^ 2 (Q k ,R k ) + A¥+V 2 (i r>k , j r>k ) + A^ +1 ' 2 {i s , k , j, tk ) 

i+l/2j + l/2 _ i-l/2.j + l/2 

+ A^ 2 (Aa k , a k ) - (5.48) 

e 

,,i+l/2,j ,,*+l/2,j , Tr i+l/2,j, \ . lTf i+l/2j/ n 



h 2 
= Al+V^fTTfc.Wfc) +A;+ 1 /2.i(Q fc) i ?fe ) +A;+V2,i (ir>fc)ir)fe) +A; +1 /2.i( is!fc) j Sifc ) 

i+l/2,j + l/2 _ 1+1/2J-1/2 

+ Aj+W(Aa ik , at ) - ^ ^ (5.49) 

Equations (4.261 - 4.266) 2 (5.50) 

i+l,j+l/2 ij + 1/2 i+l/2j" + l J+1/2J 

"* + ~ Vk = (5.51) 



6 Numerical Results 

To corroborate our theoretical results, we have implemented the variational update equations (5.39-5.43), (5.44- 
5.47), (5.48-5.51) for MHD, nematic liquid crystal flow, and microstretch continua, respectively, on a 2- 
dimcnsional cartesian grid for a variety of test cases typically used in the literature. Owing to the prevalence 
of well-established numerical test cases for MHD relative to complex fluid flow, we focus primarily on com- 
paring our novel MHD integrator with existing MHD integrators, and we give only a few examples of our 
complex fluid flow integrators for proof of concept. 

6.1 MHD Test Cases 

6.1.1 Fluid Vortex in an Azimuthal Magnetic Field 

As our first test case, we consider the evolution of a localized fluid vortex coupled with an initially azimuthal 
magnetic field in a rectangular box with vanishing boundary conditions on the normal components of v and 
B. Our setup is adapted from a classic test case from computational fluid dynamics: the advection of a vortex 
in isentropic flow [43]. We modify the setup by imposing no-transfer conditions across the boundaries of 
the domain, incompressibility of the initial velocity field, and uniformity of the fluid density. The numerical 
details of the setup are given below: 

• Domain: [0, 10] x [0, 12] 

• Boundary Conditions: Tangential Velocity Field and Magnetic Field 

• Resolution: 20 x 24 



2 Although these equations do in fact admit cartesian representations, we do not modify them for they are already in a form 
most suitable for efficient implementation as explicit updates. 
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Figure 6.1: (a) Energy vs. iteration number for an MHD simulation on a 2-dimensional cartesian grid using 
the variational integrator (5.39-5.43). (b) Cross-helicity vs. iteration number for the same MHD simulation. 
Notice that this integrator preserves the discrete cross-helicity (4.84) exactly. 



• Time Step: h = 0.5 

• Time Span: < t < 80 

• Initial Conditions: 



u[x, y) = u Q + — exp I — — 



(y - ya) 



v(x,y) 
B x (x,y) 



& n- 



"KX\ 

" 10.) m> Vl2 



(X - Xq) 



By(x,y) = cos J si 



firy 



p(x,y) 



1 - 



8771 



exp(l 



ith 



r = y/(x- x ) 2 + (y- y ) 2 
• Parameters: xq — i,ya = 5.5, uq = 0.5,/? = 5,7= 1.4. 
Since the initial velocity field defined above is inconsistent with the boundary conditions, we have used 



a numerical root-finding algorithm to obtain a nearby field that satisfies v • n 



dM 



0. 



As Fig. 6.1 shows, the integrator exhibits excellent energy behavior and preserves the discrete ana- 
logue (4.84) of the continuous cross-helicity (v b ,B) exactly. Moreover, the divergence of the velocity field 
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and the divergence of the magnetic field automatically vanish (modulo numerical roundoff errors) at each iter- 
ation. An animation of this MHD simulation may be viewed at http : / /www . its . caltech . edu/~egawlik/ 
MHD/MHDvortex . avi. 

6.1.2 Magnetic Reconnection 

A good testament to the structure-preserving nature of a numerical MHD scheme is its ability to conserve the 
topology of the magnetic field lines over the course of a simulation. Indeed, the pure advection of the magnetic 
field in ideal MHD gives rise to topological invariants that eliminate the possibility of magnetic reconnection, 
a phenomenon whereby magnetic field lines break and recombine, altering the mutual connectivity of the field 
lines. While such processes occur naturally in resistive MHD, where a magnetic resistivity term proportional 
to V 2 B appears on the right-hand side of the magnetic field evolution equation (4.55), the appearance of 
magnetic reconnection in an ideal MHD simulation is evidence of artificial numerical resistivity introduced 
by the numerical scheme. 

Here we consider a setup proposed by Gardiner and Stone [20] in their studies of numerical methods for 
compressible MHD, whereby a pair of current sheets is juxtaposed with a perpendicular velocity field on a 
periodic domain. The details of the setup are as follows: 

• Domain: [0, 2] x [0, 2] 

• Boundary Conditions: Periodic 

• Resolution: 30 x 30 

• Time Step: h = 0.1 

• Time Span: < t < 8 

• Initial Conditions: 



u(x,y) 


= u sm(ny) 


v(x,y) 


= 


B x (x,y) 


= 






1 Bo if x < X\ 


B y (x,y) 




. —Bo if xi < x < x-i 






1 Bq if X2 < x 



p(x,y) = 0.1 

• Parameters: x\ — 0.5, x^ — 1-5, uo = 0.1, Bq = 1, 9 = tan _1 (0.5). 

Snapshots of the magnetic field lines over the course of the simulation are displayed in Fig. 6.2. Notice 
that our integrator exhibits virtually no artificial reconnection, even on a relatively coarse grid. The reader 
is invited to compare with Fig. 12 of Gardiner and Stone's study [20] (copied here in Fig. 6.3), where 
the application of a constrained-transport Godunov scheme has led to visible reconnection in the form of 
magnetic field islands along the current sheets, despite their use of a high-resolution 256 x 256 grid. 



6.1 MHD Test Cases 61 



6.1.3 Field Loop Advection 

For sufficiently weak magnetic fields, the ideal MHD equations (4.54-4.57) may be viewed as a nearly decou- 
pled system whereby the magnetic field travels passively with the fluid and exerts negligible influence on the 
flow. Here we consider the advection of a weak magnetic field loop in the presence of a uniform velocity field 
on a periodic domain. The setup, proposed by Gardiner and Stone [20], is designed in such a way that the 
magnetic field loop, centered at the origin at t = 0, returns to the origin at fixed intervals of length At = 1. 
The details of our simulation are as follows: 

• Domain: [-1,1] x [-0.5,0.5] 

• Boundary Conditions: Periodic 

• Resolution: 128 x 64 

• Time Step: h = 0.01 

• Time Span: < t < 2 

• Initial Conditions: 

u(x, y) — uq cos(#) 
v(x,y) = u sin(6») 
A(x, y) = A {R - y/x 2 + y 2 ) 
p{x,y) = l 

• Parameters: v = A = 10~ 3 ,i? = 0.3,8 = tan" 1 (0.5). 

In the specifications above, A(x, y) corresponds to the magnetic vector potential, related to the magnetic 
field via B x = §|, B y = 

In Fig. 6.4, we display snapshots of the magnetic field lines at t = 0, t = 1, and t — 2. The field loop 
retains its structure and returns to the origin at the expected instants in time, although it suffers some 
distortion, presumably due to our magnetic update scheme's (5.41-5.42) low order of accuracy and absence 
of an upwind bias. Note, however, that our scheme introduces no artificial decay in the volume-averaged 
magnetic pressure B x + By during the advection of the field loop, in stark contrast with the power-law decay 
in magnetic pressure observed by Gardiner and Stone [20] in their Godunov scheme simulations of the same 
system. 

6.1.4 MHD Rotor 

As a model for star formation, the compressible MHD rotor provides an intriguing test bed for numerical 
simulation. The classic MHD rotor setup consists of a dense rotating disk of fluid in an initially uniform 
magnetic field. Rotation of the fluid induces winding of the magnetic field lines, generating magnetic tension 
that ultimately leads to a reduction in the disk's angular momentum due to the emission of torsional Alfven 
waves. Here we adapt the setup for the MHD rotor problem studied by Balsara and Spicer [6] as follows: 

• Domain: [0, 1] x [0, 1] 

• Boundary Conditions: Periodic 
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• Resolution: 30 x 30 

• Time Step: h = 0.003 

• Time Span: < t < 0.36 

• Initial Conditions: 



-uaf(r)(y - 0.5)/r if r < r 
u(x,y)={ -u o f(r)(y-0.5)/r if r < r < n 

if ri < r 

u f(r)(x - 0.5)/r if r < r 
^0~,2/) = ^ u /(r)(a;-0.5)/r if r < r < r x 

if r\ < r 



B x (x,y) = 5/V4tt 
B y (x,y) = 
p(x,y) = 1 



with 



/(r) 



VO-0.5) 2 + (y-0.5) 2 
ri — r 



• Parameters: u = 2,r = 0.1, ri = 0.115, = tan _1 (0.5). 

An animation of the evolution of the velocity and magnetic fields for this MHD rotor simulation may be 
viewed at http://www.its.caltech.edu/~egawlik/MHD/MHDrotor.avi. Snapshots of the magnetic field 
lines, displayed in Fig. 6.5, highlight the manner in which the magnetic field is "frozen" into the fluid. 

6.1.5 Orszag-Tang Vortex 

In this example, we study the evolution of current sheets in an Orszag-Tang vortex, a classic test for two- 
dimensional MHD codes. This particular example illustrates the growth of MHD current sheets: curves along 
which the magnitude of the current density j = V x B is large due to the presence of antiparallel magnetic 
field lines. Here we consider the incompressible Orszag-Tang vortex studied by Cordoba and Marliani [15]: 

• Domain: [0, 2?r] x [0, 2?r] 

• Boundary Conditions: Periodic 

• Resolution: 64 x 64 

• Time Step: h = 0.01 

• Time Span: < t < 0.75 
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• Initial Conditions: 



ip(x,y) 
A(x,y) 



2sin(y) 
cos(2y) 
1 



2 cos(x) 
2cos(a;) 



Once again, in the specifications above, A(x,y) corresponds to the magnetic vector potential, related to 
the magnetic field via B x — B y = — §^r- Similarly, ?p(x,y) corresponds to the fluid stream function, 

related to the fluid velocity field via u = , v = — ^ . 

A comparison of our integrator's performance (Fig. 6.6) with a reference solution computed by Cordoba 
and Marliani [15] (Fig. 6.7) illustrates that our integrator produces qualitatively accurate results, even at a 
low level of resolution. For comparison, Cordoba and Marliani's plots [15] were generated with an adaptive 
mesh refinement scheme using an initial resolution of 1024 x 1024; by the end of the simulation, the adaptive 
scheme employs an additional five levels of hierarchical refinement. Our simulation employed a modest 
64 x 64 grid, only lacking resolution toward the final stages of current sheet development. 

6.1.6 Convergence Rate Analysis 

As a final MHD numerical test, we study the convergence rate of our MHD variational integrator (5.39-5.43). 
We use initial conditions identical to those in Section 6.1.5, but with varying spatial resolutions (8 x 8, 16 x 16, 
32 x 32, 64 x 64, and 128 x 128) and a fixed ratio of spatial to temporal resolutions (e/h « 7.85). Fig. 6.8 
plots the L2-norm of the difference in velocity fields at t = 0.25 between successive refinements of the grid, 
as well as the difference in magnetic fields between successive refinements of the grid. The plot confirms the 
linear accuracy of the integrator. 

6.2 Complex Fluid Test Cases 

6.2.1 A Nematic Liquid Crystal Test Case: Diffusion of a Gyrating Disk 

As a numerical test of our nematic liquid crystal integrator, we consider a disk of gyrating directors immersed 
in a steady state fluid flow: 



• Domain: [0, 10] x [0, 10] 

• Boundary Conditions: Tangential velocity field, vanishing director gradient 

• Resolution: 10 x 10 

• Time Step: h = 0.4 

• Time Span: < t < 50 
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• Initial Conditions: 



. /irx\ . (iry\ 
sin — sin — 

\ 10/ V io / 



a(x,y) 
p{x,y) 



I 1 ifr<f 
[ if r > | 




1 



with 



r = 



v/(x-5) 2 + (y-5) 2 , 



where '0 is related to the velocity field (u, v) via u = v — — 

Plots of energy and angular momentum vs. time for this simulation are given in Fig. 6.9. The integrator 
exhibits good energy behavior and captures the qualitative dynamics of the system convincingly, as evidenced 
in the following animation: http : //www. its . caltech. edu/~egawlik/Complexy o 20Fluids/nematic . avi. 

6.2.2 A Microstretch Fluid Test Case: Diffusion of a Gyrating Disk 

To test our microstretch continua integrator, we consider precisely the same setup as above; the evolution 
of the system will differ from the nematic liquid crystal case due to the extra flexibility in the length of the 
director. The setup is generalized as follows: 



• Domain: [0, 10] x [0, 10] 

• Boundary Conditions: Tangential velocity field, vanishing director gradient 

• Resolution: 10 x 10 

• Time Step: h = 0.4 

• Time Span: < t < 50 

• Initial Conditions: 




R(x,y) 
a(x,y) 
X(x,y) 
jr(x,y) 
js{x,y) 

p{x,y) 



o 



















1 



with 



v / (x _5 ) 2 + (y _ 5)2; 
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where tp is related to the velocity field (u, v) via u = ^ , v = — ^j. 

Plots of energy and angular momentum vs. time for this simulation are given in Fig. 6.10. Like the 
previous integrator, the microstretch fluid flow integrator exhibits good energy behavior and captures the 
qualitative dynamics of the system convincingly; see http://www.its.caltech.edu/~egawlik/Complex7o 
20Fluids/microstretch. avi for an animation. Notice the heightened intricacy in the dynamics of mi- 
crostretch continua over nematic liquid crystals owing to the flexibility of the director field. 

7 Conclusions and Future Directions 

This study has merged techniques from variational mechanics and discrete exterior calculus to design struc- 
tured integration algorithms for magnetohydrodynamics and complex fluid flow. These integrators exhibit 
several desirable features that distinguish them from conventional integration schemes: they are symplectic, 
exhibit good long-term energy behavior, and conserve momenta arising from symmetries exactly. Moreover, 
their formulation in the framework of DEC ensures that Stoke's theorem holds at the discrete level, leading 
to an automatic satisfaction of divergence-free constraints on the velocity and magnetic fields in the resulting 
numerical schemes. 

Numerical tests of these structured integrators on cartesian grids confirmed our theoretical results and 
indicated that in addition to the properties listed above, the integrators faithfully respect topological features 
of the flow in MHD (e.g. the forbiddance of magnetic reconnection) and are predictive at remarkably coarse 
resolutions, despite being of low order. 

In the process of designing structured integrators for MHD and complex fluid flow, we uncovered a gen- 
eral framework for the design of variational integrators for continuum systems described by Euler-Poincare 
equations with advected parameters, a particular class of differential equations obtained from a variational 
principle. This finding is tremendously encouraging given the plethora of physical systems arising in mechan- 
ics, fluid dynamics, and plasma physics that are governed by such equations; see references [27] and [35] for 
examples. To name one such example, nearly all of the governing equations for geophysical fluid dynamics — a 
field that plays an instrumental role in the modeling of global climate change — are of this type [26]. We en- 
vision the exterior calculus based variational framework serving as a cornerstone for the design of structured 
integrators for these families of physical systems. 

Future Work. A number of numerical and theoretical challenges are left unaddressed in this study. 
Incorporation of compressibility, viscosity, resistivity, and inhomogeneity into the discrete variational formu- 
lation of fluid flows are obvious extensions. We are also developing a more general approach to handle Lie 
group valued functions in arbitrary dimension to resolve the issue mentioned in Section 4.2.3. On the prac- 
tical side, we have yet to develop high order versions of our variational update equations valid on arbitrary 
grids. Finally, we are interested in developing other structure- preserving integrators for, e.g., relativistic 
physical systems. 
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A Appendix 

A.l Lie Groups and Lie Algebras 

In this section, we present those aspects of Lie group theory most relevant to our development of discretiza- 
tions of continuum theories. The notations appearing here are used consistently throughout the paper. For 
a more complete treatment of Lie groups, Lie algebras, and group actions, the reader is referred to [1]. 

Let G be a Lie group; that is, G is a smooth manifold equipped with a group structure for which the 
group multiplication and inversion are smooth maps. The Lie algebra g of G is by definition the tangent 
space of G at the identity e G G: 

g=T. G ={eie=| 



g(i) for some curve g : (-6, b) C K -> G with g(0) — e f (A.l) 

t=o 



We denote the space dual to g by g*; this is the space of linear maps from g into K. We use angular brackets 
to denote the natural pairing (■,■): g* x g— > R. 

We use the following notation for left and right translation of tangent vectors in TG: 



TL gV = t_ 



TR gV = 



gh{s) (A.2) 
h(s)g, (A.3) 



where g 6 G, n £ TJ^G, and h : (—6, b) C K — > G is a curve in G with h'(0) = rj. If there is no danger of 
confusion, we simply denote these maps via concatenation: gr\ — TL g rj and r\g — TR g rj. 

A. 1.1 Group Actions 

Let V be a vector space. A left action or left representation of G on V is a smooth map <!> : G x V — > V for 
which 

1. $(e, a) = a for every a E V. 

2. $(5, a)) = $(5/1, a) for every g,h e G,a e V. 

We denote such an action via concatenation, or if there is danger of confusion, we insert a binary operator: 

ga = $(<?, a)=g-a. (A.4) 
The induced infinitesimal action of g on is given by differentiating the action of G on V at the identity: 

d 



e ' a= dt 



g{t) ■ a, (A.5) 



where g : (—6, b) C K — > G is a curve in G with g(0) = e, g'(0) = £. 

A corresponding notion is that of a right action or right representation of G on V, which is a smooth 
map \& : V X G — > V for which 

1. ^(a, e) = a for every a£ 

2. *(*(a, g),h) = *(a, for every g,h <E G,a EV. 
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We denote such an action via concatenation or with a binary operator: 

ag = ®{a,g) = a ■ g. (A.6) 
The infinitesimal right action of g on V is given by 

a-g{t), (A.7) 

t=o 

where g : (—6, b) C M. — > G is a curve in G with g(0) — e, g'(0) = £. 

All of these notions make sense if V is replaced by a manifold M; in such a scenario, we say that there 
is an action of G on M. (The term representation, however, is reserved for group actions on vector spaces 
only.) 




A. 1.2 Adjoint and Coadjoint Actions 



For any Lie group G, there is a natural representation Ad : G x g — > g of G on its Lie algebra g called the 
adjoint action, given by conjugation: 



Ad gV = TL„ ■ TR q -i ■ r] 



d 
ds 



(A.8) 



s=0 



where g € G and ft- : (~b, b) C M — > G is a curve in G with /i(0) = e, h'(0) = -q. Its dual is denoted Ad* and 
is referred to as the coadjoint action of G on g* when g is replaced by g~ x . Explicitly, 



(Ad>,r/) = {n,A.d g rf) 



(A.9) 



and the coadjoint action for fixed g £ G is the map Ad*-i : g* — > g*. The inversion of g is a convention 
ensuring that this action is a left action on g* rather than a right action. 
The bracket, or commutator, on g is the differential of the adjoint action: 



dt 



Ad 



t=o 



dt ds 



g(t)h(s)g(ty 



(A.10) 



s,t=0 



where g and ft are curves in G passing through the identity at t = and with </(0) = ^, ft'(O) = r\. 
The infinitesimal adjoint action is the map ad : g x g — > g defined by 

ad € r? := %r}\. 



(All) 



This notation permits the definition of the infinitesimal coadjoint action ad* of g on g* , given by the dual 
of the infinitesimal adjoint action: 

(ad^,?7) = {/jL,ad^}. (A.12) 

A. 1.3 The Exponential Map 

Given any Lie group G, there is a unique map exp : g — > G satisfying 



| exp(tO = TR cxm) tt = TL cxpm t£ 
exp(O) = e 



(A.13) 
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for all £ € g, t G R. Wc call this map the exponential map. It is shown in standard texts on differential 
geometry that the exponential map is a local diffcomorphism that maps a neighborhood of £ g to a 
neighborhood of e e G. 

For matrix groups, the exponential coincides with the usual matrix exponential: 

(tA) 



exp(M)=£^. (A.14) 



^ 

z=0 



A. 2 Exponentiation on the Nematic Liquid Crystal Configuration Space 

This section presents proofs of four lemmas stated in Section 4.2.3, where a study of exponentiation on the 
spatially discretized nematic liquid crystal configuration space G = T>(M) (S) ^°(M, S 1 ) was conducted. We 
restate the lemmas themselves for the reader's convenience. 

Lemma, (cf. Lemma 4.1) The exponential map exp : g — > G for the group G = D(M) (S) il®(M, S 1 ) is given 
by 

exp(t(A,w)) = (e^^A-^I -e- tA )io), (4.142) 

where t £ K 7 A € D(M), ui G f2°(M), etA * s ^ e usual matrix exponential, and A _1 (/ — e~ tA ) is to be regarded 
as a power series 

A- 1 (I-e- tA )=tI- 



t 2 A t 3 A 2 



2 6 

(so that it is defined even for A not invertible). 
Proof. Set t — to check that 

exp(0,0) = (7,0). 
Now differentiate (4.142) with respect to t to confirm that 

4- expMA, u)) = t(Ae tA , e - tA u) 
at 

= TR C x P (t(A,u)){tA,tu}) 



Lemma, (cf. Lemma 4.2) The map r : g —¥ G given by 

t(A,u) = (cay(A),cay(-A/2)w) (4.144) 

is a group difference map for the group G = P(M) (§)f22(M, S 1 )- That is, t is a local approximant to exp 
with t(0, 0) = (I, 0), and r satisfies 

t^lo)- 1 =t(-A,-w) V(4,w)ej. (4.145) 

Proof. Comparison of the series expansion of (4.142) with that of (4.144) confirms that r is a local approx- 
imant to exp. The proof of (4.145) is by direct calculation. ■ 

Lemma, (cf. Lemma 4.3) The inverse right-trivialized tangent dr^ 1 : g x g — > g of the map (4-144) * s given 
by 

dT {A,u)( B ' VO = ( dcay^B, cay(-A/2)V + -<k&y A/2 (dc&y-^B)w j . (4.146) 
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Proof. By formulae (3.28), (4.142), and (4.144), 
d 



de 
d_ 
de 
d 
de 
d 
de 



t~ (exp(s(B, 4>))r(A, oj)) 



— 



- 1 {{e sB ,B-\l-e- EB )i>){ C zv{A), cay(-A/2) W )) 



e=0 



t- 1 {e sB cay{A), cay(- A)B~ 1 (I - e~ sB )^ + cay(-A/2)w) 



6 = 



(cay- 1 (e £B cay(A)), 



6 = 



cay 



cay _1 (e e -°cay(A)) 



|cay(-A)S _1 (/ - e~ sB )i/j + cay(-A/2)w 



dcay/B, cay(A/2)[cay(-A)^] + -dca,y A/2 (dca,y A 1 B)uj 



dcay/B, cay(-A/2)^ + -dcay j4/2 ( ( icay A 1 B)t 



Lemma, (cf. Lemma 4.4) For the map r given by (4-144), the dual of the operator dr 1 is given by 

(^(A^n^) = (Wy^ 1 )*^ + ^(dcay^ 1 )*(rfcay A/2 )*skew(7r W T ), cay(A/2)^ . 
iVoie i/iai if t: is a scalar multiple of lo, then this reduces to 

^(L))*^^) = (W^ 1 )*^, cay(A/2)7r) . (4.148) 
Proof. One checks that the above formula satisfies 



A. 3 Derivations of Cartesian Realizations 

This section presents derivations of the results reported in Section 5. We employ the notation of Section 5 
and make use of Fig. 5.1 throughout the discussion. 

A. 3.1 The Lie Derivative of a Discrete One- Form 

Let A, C £ S be discrete vector fields that approximate continuous vector fields u = (it, t>),w = (r, s) £ 
%div(M), respectively, as in Section 5.1.1. 

The entries of the matrix C b G fl^M) / dfl° d (M) are then related to those of C via (2.29). Note the 
disparity between the sparsity patterns of C and : C has nonzero entries associated with cell-to-neighbor 
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transfers only; C b has nonzero entries for cell-to-neighbor transfers and for transfers between cells and their 
neighbors' neighbors. 

Fix a pair of vertically adjacent cells C m and C n . For I € {1, 2, . . . , N}, let N(l) denote the set of cells 
sharing an edge with the cell Ci . Taking into account the sparsity patterns of A and C b , we have 

(■£^(7 ) mn = [C\^]mn 

= ^ C b nl Ai n — ^2 A m \C\ n . (A. 15) 

lGN(N(m))r\N{n) l£N(m)nN(N(n)) 

Now, with the aid of Fig. 5.1(b), substitute the values of the entries of C and A in accordance with (5.10-5.13) 
and (2.29): 

(£ A C') mn =-\{ s i+1 ^ + j+W+A v i+1 ^ +1 



i+\/2,j-l , „i+l/%j\ „i+l/2,3-l 



2 
1 

2 

1 / eu l+l ' 3 (r, S) i+1 j + 1 /2 i+l/2j \ i+lj + l/2 



1 / eu l+l > 3 (r, s) _ i+lij _ 1/2 l+ i/ 2 ,j\ i+ij-x/2 
1 / ew^{r, 3 ) i;i+1/2 _ i+1/2 ,A i)j+ i /2 



1 ^(r,,) _ riJ _ 1/2 _ s i+i/ 2 ,A u „_ 1/a> : A .l«i:, 



2 

Here, we have expressed the entries of C in terms of the quantity u>(r, s) in anticipation of the role played 
by the discrete curl of the vector field w = (r, s): 

r M-l/2 , s i+l/2,j _ r i,j+l/2 _ s i-l/2,3 

ui l ' J (r,s) = (A.17) 

Rearranging, we obtain 

(^CVn = - 2 K J s ) ( 2 J + w J ( r > s ) ( — 2 - ) 



_ I s i+l/2,j ^i+l/2,j+l _ w 'i+l/2j + u i+lj + l/2 _ u ij + l/2^j 
_ I s i+1/2J ^i+1/2,3 _ w i+l/2,j-l + u i+l,i-l/2 _ u i,i-l/2^ 

_ I ^ sv y+i/2, j+ i + ^y+id+i/2 + ( 8V y+i/2J + (ruy^+y^ 

+ I ((to^+Vsj + (ru )i+W-i/2 + ^m/2,,-1 + {ru y- 3 -i/2^ (A lg) 

The second and third terms in the expression above are zero by the discrete divergence-freeness of u. The 
last two terms constitute a full discrete differential. Thus, modulo a discrete differential, 

{£ A &) mn = -- i^(r,s) I j +^ +1 ^(r,s) f j j (5.14) 
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Proceeding similarly for a pair of horizontally adjacent cells C m and C n centered at (i — 1/2, j + 1/2) and 
(i + 1/2, j + 1/2), one obtains 



(£ A C b ) 



V *-1/2J _|_ w i+l/2,j 



(r, s) 



+ ^+1/2,^+1 



(5.15) 



A. 3. 2 The Lie Derivative of a Discrete Vector Field 

Let A, B £ S be discrete vector fields that approximate continuous vector fields u = (u,v),v = (p,q) £ 
Xdi v (M), respectively, as in Section 5.1.2. Below we derive the cartesian realization of £ aB = —[A,B] £ 
[S,S\. ' 

Before doing so, recall from Section 2.3 that weak equalities involving elements of [S,S] have solutions 
in S that are expressible in terms of the sparsity operator * defined in (2.35). Let us therefore compute the 
entries of the matrix £aB — — [A,B] £ [S,S] and apply the sparsity operator (2.36) to recover a matrix 
(£aB)^ £ S that is weakly equal to £aB. 

Note that for any cells C a and Cf, separated by a distance of two, 



(£AB)ab 



E 



B a cA c b — A ac B c 



(A.19) 



c£N(a)nN{b) 



Thus, in terms of the vector fields u = (it, v) and v = (p, q) corresponding to A and B, respectively, the (m, n) 
entry of the matrix (£aB)^ for any pair of vertically adjacent cells C m and C n centered at (i + 1/2, j — 1/2) 
and (i + 1/2, j + 1/2) is given by 



(£AB) mn — 



4e 2 



q i+1/2 ' j (u iJ + 1/2 + V i+1/2J - U l+1 " 3 + 1/2 ) - t, i+1 /2,j( 2 /J + l/2 + q i+l/2,j _ p i+l,j + l/2^ 
+ (q i+1/2 ' j +p i+1 'J- 1 /2 _ p hi-l/2\ v i+l/2,j _ + ^i+lj-1/2 _ ^ij-l/as i+1/2 ,j 

+ I{gi+l/2j" u i+lJ+l/2 _ v i+l/2,j p i+l,j+l/2 +p i+~L,j-l/2 v i+3/2,j _ ^+1^-1/2^+3/2^ 
+ p i,j-l/2 v i+l/2,j _ B i,j-l/2 i+l/2,j + i-l/2,j u i,j+l/2 _ ,j i ,3+1/2 

_|_ i+lj + l/2 w i+3/2,j _ M +lj+l/2 i+3/2j + t+l/2,i u i+X,j-l/2 _ ^i+l/2,j -i+lj-1/2 



+p hj + l/2 v i+l/2,j _ u i,j+l/2 q i+l/2,j + q i-l/2,j u i,j-l/2 _ j p i,j-l/2\ 



(A.20) 



Note that in the first two lines of the equation above, we invoked the discrete divergence-freeness of 
u = (u, v) and v = (p, q) to rewrite the fluxes across the top center and bottom center edges in Fig. 5.1(b) 
in terms of the fluxes across labeled edges in the figure. 

Upon simplification, (A.20) reduces to 



(£aB)\ 1ui — 



1 f 1 



2e e 



u i+l,j-l/2 + u i+l,j+l/2 

2 

u i,j-l/2 +u i,3+l/2 
2~ 

i+lj-1/2 +p i+l,j+l/1 
2 



Jl+3/2,3 + qi+1/2,3 

2 

q i+l/2,j + q i-l/2,j 

2~ 

v i+3/2,j + v i+X/2,j 

2 



M-l/2 +p i,j + l/2\ / v i+l/2,j +v i-\/2,j 



(A.21) 
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In the abbreviated notation discussed in this Section 5's introduction, (A. 21) reads 

{£A Bf mn = -1 ^^-W* _ p^yW-W^ (5 21) 

Proceeding similarly for a pair of horizontally adj acent cells C m and C n centered at (i — 1/2, j + 1/2) and 
(i + 1/2, j + 1/2), one obtains 

i / ^j+yj+i _ M i, Jg «,j g+yj+i-j^ \ 

A. 3. 3 The Antisymmetrization of an Outer Product 

Let E,F € Oj(M) be discrete zero-forms approximate continuous scalar fields a, j3 G F(M), respectively, as 
in Section 5.1.3. 

Let us compute the (m, n) entry of skew(F£' T ) for vertically adjacent cells C m and C n centered at 
(i + 1/2,. j - 1/2) and (i + 1/2, j + 1/2), respectively: 

skew(F£ T ) mn = l - ({FE T ) mn - (EF T ) mn ) 

= i( ( a i + 1 /2,i-l/2 Q ,i+l/2,j+l/2 _ a t+l/2,j-l/2 j gt+l/2,:7+l/2N 

^+1/2,^-1/2 + Q i+l/2 J + l/2 \ / pi+l/2,j+l/2 _ ^1+1/2,^-1/2 

"1 / l~ ^ 

^i+l/2,j-l/2 + ^i+l/2j + l/2\ / a i+l/2j + l/2 _ Q ,»+l/2j-l/2 

2 

, / . /flt+l/2,i+l/2 _ m+l/2,j-l/2\ / i+l/2j+l/2 _ i+l/2,j-l/2 
e ' '-^1/2,3 / P P \ _ oi+l/2,j I « « 



= -2 [ a t+i " 3 [ - - ) - a 



(5.27) 



Proceeding similarly for a pair of horizontally adjacent cells C m and C n centered at (i— 1/2, j + 1/2) and 
(?' + 1/2, j + 1/2), one obtains 

e/ , / fl*+l/2,J+l/2 _ gi-l/2,j+l/2 \ / i+l/2,j+l/2 _ ^-1/2,3+1/2 

skew^U = -| U«+V= ^ -P j - (<■ 

(5.28) 
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Figure 6.2: Evolution of the magnetic field lines for an MHD current sheet simulation on a 2-dimcnsional 
cartesian grid using the variational integrator (5.39-5.43). The frames above are snapshots at times t = 
0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, in normal reading order from left to right, top to bottom. 
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Figure 6.3: (Borrowed from Gardiner and Stone [20].) Evolution of the magnetic field lines for an MHD 
current sheet simulation on a 2-dimensional cartesian grid using a constrained transport Godunov scheme. 
The frames above are snapshots at times t — 0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0, in normal reading order 
from left to right, top to bottom. Compared with our variational integrator (see Fig. 6.2), the Godunov 
scheme exhibits a noticeably larger amount of artificial magnetic reconnection due to numerical resistivity. 
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Figure 6.4: (a) Evolution of the magnetic field lines for an advected magnetic field loop using the variational 
integrator (5.39-5.43). The frames above, from top to bottom, are snapshots at times t = 0, 1, 2. (b) Decay 
of the volume-averaged magnetic pressure B% + By during advection of the field loop. Unlike the Godunov 
scheme proposed by Gardiner & Stone [20], our integrator exhibits no noticeable decay in net magnetic 
pressure during what amounts to a virtually passive advection of the magnetic field loop. 
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Figure 6.5: Evolution of the magnetic field lines for a simulation of the incompressible MHD rotor on a 
2-dimcnsional cartesian grid using the variational integrator (5.39-5.43). The frames above are snapshots at 
times t = 0.000, 0.042, 0.084, 0.126, 0.168, 0.210, 0.252, 0.294, 0.336, in normal reading order from left to right, 
top to bottom. 
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Figure 6.6: Current density contours for a simulation of the incompressible Orszag-Tang vortex on a 2- 
dimensional cartesian grid using the variational integrator (5.39-5.43). The frames above are snapshots at 
times t = 0.00,0.25,0.50,0.75, in normal reading order from left to right, top to bottom. 
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Figure 6.7: (Borrowed from Cordoba and Marliani [15].) Current density contours for a simulation of the 
incompressible Orszag-Tang vortex on a 2-dimensional cartesian grid using an adaptive mesh refinement 
scheme. The frames above are snapshots at times t = 0.00, 0.25, 0.50, 0.75, in normal reading order from left 
to right, top to bottom. 
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Figure 6.8: Convergence rates for a simulation of the incompressible Orszag-Tang vortex on a 2-dimensional 
cartesian grid using the variational integrator (5.39-5.43). Data points represent the L2-norm of the difference 
in vector fields at t = 0.25 between successive power-of-two refinements of the grid. The reported grid 
resolution at each data point corresponds to the finer of the two resolutions under scrutiny. 
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Figure 6.9: (a) Energy vs. iteration number for a nematic liquid crystal simulation on a 2-dimensional 
cartesian grid using the variational integrator (5.44-5.47). (b) Total angular momentum vs. iteration number 
for the same nematic liquid crystal fluid simulation. Notice that this integrator preserves the total angular 
momentum exactly. 
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Figure 6.10: (a) Energy vs. iteration number for a microstretch fluid simulation on a 2-dimensional cartesian 
grid using the variational integrator (5.48-5.51). (b) Total angular momentum vs. iteration number for the 
same microstretch fluid simulation. Notice that this integrator preserves the total angular momentum exactly. 



